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mini-j-  of  any  twice  continuously  differentiable  real  function  of  n  real 
variables  on  a  closed,  bounded  domain  is  described.  The  algorithm  also  provides 
infallible  bounds  on  the  location  of  the  global  minimum.  The  aloorithr  uses 
Interval  arithmetic  and  requires  the  availability  of  several  fundamental 
Interval  arithmetic  processors  for  its  operation. 


toiTio**  o»  >  ncv  *»  >*  oatoi.CTc 


UNCLASSIFIED 


WlOBALM^N  is  a  program  for  finding  the  global  minimum  of  a  nonlinear 
function  f  C  of  n  variables.  It  provides  infallible  bounds  on  the  minimum 
value  F*  of  f  in  any  prescribed  box  (  a  box  is  a  parallelepiped  with  sides 
para'Ie'  to  the  coordinate  axes).  It  also  provides  infallible  bounds  on  the 
point  s}  x*  at  which  the  global  minimur  occurs.  t  t1  ) 

It  is  assuned  that  the  minimum  occurs  at  a  stationary  point  of  f.  It 
is  possible  to  modify  the  program  to  allow  for  the  case  in  which  the  global 
mini-jr  is  not  a  stationary  point  and  hence  occurs  on  the  boundary.  The 
necessary  modifications  are  described  in  Appendix  A,  which  is  to  appear  in 
Nurerische  Vathematik.  This  paper  describes  the  theoretical  and  practical 
considerations  used  in  writing  GL0BALK1N.  It  also  provides  the  results  of 
severa'  test  problems  treated  by  GLOBALMIN. 

. ne  work  described  here  directly  builds  upon  earlier  work  carried  out 
at  Washington  State  University  in  the  general  area  of  interval  arithmetic. 

In  particular,  three  reports  by  Carol  Clarke  that  describe  fundamental  pro¬ 
cessors  used  by  Gl06ALMI*»  are  included  as  Appendices  ir.  this  report  for  com¬ 
pleteness.  We  describe  particular  considerations  in  regard  to  this  earlier 
work  that  are  germane  to  GLOBALMIN  in  the  next  two  sections. 
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2.  THE  CIUDGE  PRECOMPILER 


GLQBALMIN  was  written  assuming  a  precompiler  named  CLUDGE  would  be  used 
The  main  purpose  of  CLUDGE  is  to  accept  and  work  with  intervals  in  the  same 
way  Fortran  works  with  real  numbers,  complex  numbers,  integers,  etc.  It 
also  provides  interval  subroutines  for  evaluating  the  elementary  transcen¬ 
dental  functions  with  interval  arguments. 

„e  now  briefly  describe  some  of  the  features  of  CLUDGE.  For  a  more 
detailed  discussion,  see  Appendix  I. 


2.1  Interval  Variables 

In  order  for  CLUDGE  to  accept  a  variable,  say  A,  as  an  interval,  the 
variable  must  be  so  declared  as  INTRVAL  A.  It  is  possible  to  dimension  an 
interval  variable.  For  example,  SNTRVAL  A  5,7).  To  read  in  or  print  out 

O'" 

an  interval.  jO*  must  leave  roor  for  twc  nurbers.  For  example. 


101  FORMAT  (2r 10.0) 


When  you  do  use  READ  or  WRITE  statements  that  involve  Interval 
variables,  a  warning  appears.  St's  harmless  and  usually  unnecessary 


2.2  Miscellaneous  Warnings 

Although  the  CLUDGE  package  includes  many  functions,  it  doesn't  declare 
their  type  for  you.  For  example,  if  you  wish  to  use  1NTU0N,  an  interval¬ 
valued  function,  you  must  declare  INTUON  as  an  interval  name.  For  example: 


interval  intuon,  ta.tb.tc 


TC-INTUON(Ta’.TB) 

As  another  example,  INTINF  and  INTSUP  are  real-valued  functions  and 
must  be  declared  so: 

REAL  INTINF, INTSUP 
INTRVAL  TA. 

P* INTINF (TA) 

Q-INTSUPNTA) 

Finally,  when  you  write  your  own  functions,  explicitly  state  the  type, 
ever  if  it  isn't  INTRVAL. 


3.  NOTATION  AND  TERMINOLOGY 

In  the  comments  within  the  program  and  within  the  following  writeup,  the 
authors  have  used  and/or  invented  some  unusua:  terns.  A  few  comments  are  in 
order.  As  an  example,  consider  AIA(I-N.K). 

If  you  encounter  a  symbol  such  as  the  above  in  a  comment,  it  refer $  to 
the  vector  residing  in  the  consecutive  memory  locations  AIA’l.K).  AIA(2,h}... 
AIA;n,K).  for  this  particular  example,  it  is  an  interva’  vector. 

'he  comments  nak.e  use  o*  the  following  terminology: 

A  box  is  the  cartesian  product  of  N  intervals.  Equivalently,  a  box 
is  ar,  interval  vector.  Note  that  AIA(I-N.k)  is  the  h-th  box  ir,  the  list  of 
boxes  AIA. 

The  width  of  a  bo>  is  the  width  of  the  widest  interval  defining  the  box. 
Thus  if  X  is  an  interval  vector  box)  with  components  X,.  (i  *  .  .  ,N)  and 

if  w  x.)  denotes  the  width  of  the  interval  X ^ ,  then  the  width  of  the  box  is 
max  w(X^) 

1<1<N 

The  functional  width  of  a  box  is  the  width  of  the  interval  obtained  by 
evaluating  the  object  function  f  over  the  box  (using  interval  arithmetic). 
Note  that  this  interval  contains  the  range  of  the  function  f  over  the  box. 

The  driver  MAIN.  The  driver  for  the  program  is  called  MAIN.  It  is 
used  to  call  GLOBE  which  is  the  central  subroutine  of  the  program. 

The  subroutine  GLOBE.  The  minimization  algorithm  uses  a  number  of  sub¬ 
algorithms  to  find  the  global  minimum.  These  subroutines  are  called  by 

globe. 

Explanation  of  arguments.  The  calling  statement  for  GLOBE  is 


CALL  GLOBE  (N,  XL,  XR,  AIK,  AIWM,  ARK,  ARWM,  ARWMA,  AZW,  AIA, 
AIB,  AIC,  NDIM,  KENDC,  AID,  EA,  EB,  P,  FLOW,  FBAR,  NSTEP) 


The  arguments  are: 

N  is  the  dimension  of  the  domain. 

XL  is  the  vector  of  left  endpoints  of  the  current  box.  Thus  XL ( I ) 

(I  1  1 , . . . ,N)  is  the  left  endpoint  of  the  interval  defining  the  I-th  component 
of  the  box. 

XR  is  the  vector  of  right  endpoints  of  the  current  box. 

AIW  is  an  interval  matrix  of  dimension  N  by  8.  It  is  used  as  a  set  of 
interval  work  vectors. 

AIWM  is  an  N  by  N  matrix  containing  the  interval  Hessian. 

A**  is  a  real  N  by  8  matrix.  It  is  used  as  a  set  o*  real  wc  *  vectors. 

AR tF  is  a  real  N  by  N  .natrix  containing  the  midpoint  of  AIWM. 

ARa'MA  is  a  real  N  by  V  matrix  containing  the  approximate  inverse  of 

ARWK. 

AZ*  is  an  integer  vector  of  o**der  N.  It  is  used  to  store  flags  for 
various  Purposes. 

AIA  is  an  interval  matrix  with  N  rows.  It  contains  the  insufficiently 
small  boxes.  Each  column  is  a  box  generated  by  the  program  and  the  number  of 
boxes  changes  as  the  progra-  rjns.  for  small,  easy  problems  the  number  of 
colu-ns  required  may  be  1C  or  less.  For  larger  or  more  difficult  problems 
it  is  safer  to  allow  for  100  boxes.  The  program  is  written  so  that  the 
current  box,  which  is  chosen  from  AIA,  tends  to  be  the  largest  box  in  AIA. 

This  improves  the  speed  of  the  algorithm  but  causes  the  number  of  boxes  to 
be  larger.  To  reduce  storage  (while  increasing  the  run  time  somewhat), 
the  program  can  be  modified  so  that  the  current  box  is  chosen  to  be  the 
smallest  one  in  AIA. 

A I B  is  an  interval  matrix  with  N  rows.  It  contains  boxes  small 
enough  to  satisfy  the  error  criterion.  As  for  AIA,  the  number  of  columns 
depends  on  the  problem.  If  the  problem  has  only  one  or  two  solutions,  AIB 
should  never  contain  more  than  5  or  10  boxes.  If  the  problem  has  m 
solutions,  then  AIB  will  probably  contain  at  least  m  boxes  eventually. 

AIC  is  an  interval  matrix  with  N  rows.  It  contains  boxes  which  satisfy 
both  the  requirement  of  smallness  (width  EA)  and  the  criterion  of  sms'll 
functional  width  (  EB).  At  the  end  of  the  run,  all  remaining  boxes  will  be 
in  AIC.  Since  tentative  solution  boxes  are  also  put  temporarily  in  AIC,  it 


must  contain,  say,  ten  more  columns  than  there  are  solution  points.  The  pro¬ 
gram  could  be  rewritten  to  occasionally  purge  AIC  of  boxes  which  do  not  contain 
solutions.  Thus,  the  dimension  of  AIC  could  be  smaller  than  is  currently 
necessary. 

NDIM  is  the  number  of  columns  in  AIA.  Since  the  number  of  boxes  to  be  stored 
in  AIA  is  not  known  in  advance,  it  must  be  guessed.  A  safe  choice  (generally) 
is  100. 

KENDC  is  the  number  of  boxes  in  AIC. 

AID  is  an  interval  vector  containing  the  functional  widths  of  the  boxes  in 
AIC.  Thus  it  has  the  sane  number  of  elements  as  the  number  of  columns  in  AIC. 

EA  is  ar,  error  tolerance.  Any  final  box  must  be  less  than  £A  in  width. 

E5  is  an  error  tolerance.  Tne  functional  widtn  o#  f  over  each  final  box 

must  be  less  than  EB.  When  tne  algorithm  terminates,  we  will  have  FBAR-FlOW<EB. 

P  is  a  parameter  which  has  beer  set  (rather  arbitrarily)  to  0.75.  If 
during  one  complete  iteration,  the  current  bos  is  not  reduced  in  width  by  at 

least  the  factor  1-f,  then  the  box  is  split  in  half.  More  experience  is 

needed  to  choose  P  optimally. 

PLOW  is  a  lower  bound  on  the  global  minimur  f*. 

FBAR  is  an  upper  bound  on  f*. 

NSTEP  is  the  njr.ber  of  iterative  steps  used  by  algorithm.. 

User  supplied  subroutines.  The  user  must  supply  subroutines  to  evaluate 
the  object  function  f,  its  gradient  g,  and  its  hessian  H.  The  program  requires 
interval  forms  of  these  functions.  Tne  sharpness  of  the  results,  and  hence  the 
rate  of  convergence  of  the  algorithm,  depends  on  the  steps  used  to  evaluate 
these  functions.  Hence  the  subroutines  should  be  written  so  as  to  obtain  as 
sharp  results  as  feasible. 

The  object  function  f(x)  must  be  evaluated  as  an  interval  function  named 
FI  and  be  of  the  form, 

INTRVAL  FUNCTION  FI  (N.AIW.KNEXTA.L) 

Here  N  is  the  dimension  of  the  problem,  AIW  is  matrix  wherein  the  interval 
'-argument  vector  (the  box  over  which  f  Is  being  evaluated)  is  stored  and  KNEXTA 
Is  the  specific  column  of  AIW  in  which  it  is  stored.  L  is  a  flag  which  allows 
recovery  or  termination  if  difficulty  (such  as  overflow)  occurs  in  the  evalua¬ 
tion.  Set  L  *  0  if  all  goes  well.  Set  L  *  1  otherwise.  No  actual  use  of  the 
flag  is  made  by  the  program.  Its  use  is  Included  to  aid  the  user  in  recovery 
or  termination  if  he  chooses  to  modify  the  program  appropriately. 


The  interval  gradient  must  be  called  GI  and  be  evaluated  using 
SUBROUTINE  GI  (N.AIW, KNEXTA. AJW, JCOL ,L) . 

The  arguments  N,  Alw,  KNEXTA  and  t  are  the  same  as  for  FI.  The  argument 
AJW  specifies  the  matrix  in  whicn  the  interval  gradient  is  to  be  storec  and 
JCOL  specifies  the  column  of  AJW. 

"he  interval  Hessian  must  be  called  HI  and  be  evaluated  using 
SUBROUTINE  HI  (N.AIW, KNEXTA, AIWM.L). 

The  arguments  N,  AIW,  KNEXTA,  and  L  are  the  same  as  for  FI  or  GI.  The 
interval  Hessian  should  be  placed  in  the  matrix  AIWM. 

When  computing  the  Hessian,  different  input  arguments  for  different 
elements  car  be  real  degenerate  intervals).  See  the  accompanying  paper  by 
t.  Hansen.  Real  arguments  should  be  used  when  possible  to  reduce  the  widths 
of  the  interval  elements  of  the  Hessian. 

"here  is  a  need  to  have  the  diagonal  elements  of  the  Hessian  with  all 
interval  elements.  This  is  used  for  a  non-convexity  check.  Thus  the 
user  must  Supply 

SUBROUTINE  HO I  AG  (N, AIW, KNEXTA, JCOL). 

The  argument  \,  AIW,  and  KNEXTA  are  as  in  HI.  The  argument  JCOL  specifies' 
the  column  o *  AIW  in  which  the  diagonal  of  the  Hessian  (with  non-degenerate 
input]  should  be  stored  as  an  interval  vector. 

Miscellaneous  consents.  The  mode  used  does  not  allow  passing  of 
para~eters.  This  profcle"  is  circumvented  by  placing  parameters  in  a  COMMON 
block  and  referencing  the"  in  both  the  main  program  and  the  five  subprogram . 
However,  space  should  be  reserved  in  the  beginning  of  the  COMMON  block  for 
the  integers  I  FI ,  IGI,  and  IHI.  If,  in  the  main  program,  each  of  these 
counters  is  initialized  to  zero,  and,  in  each  subprogram  these  counters  are 
incremented,  one  can  easily  keep  track  of  how  many  times  each  subprogram 
has  been  called. 

Subroutines.  The  central  subroutine  in  GLOBALMIN  is  GLOBE.  It  calls 
*  various  other  subroutines  which  we  now  describe  briefly. 

JOIN  forms  an  interval  number  from  two  real  numbers. 

MCHK  determines  whether  a  given  interval  contains  zero. 


SQUARE,  SQURT,  and  CUBE  find  the  square,  square  root,  and  cube, 
respectively,  of  an  interval  number. 

CHk  examines  a  box  to  see  if  it  is  less  than  EA  in  width.  Any  such  box 
in  AIA  is  moved  to  AJB  unless  it  is  degenerate  (a  single  point).  In  the 
latter  case,  the  box  is  movea  to  AIC. 

SPLIT  divides  a  box  into  two  pieces  and  places  each  in  AIA. 

MONO  examines  the  interval  gradient.  If  a  component  is  positive  or 
negative  over  a  box,  then  f(x)  is  monotonic  in  the  corresponding  direction  and 
cannot  have  a  stationary  point  in  the  box.  Therefore,  the  box  is  eliminated. 

COEf  'inds  the  coefficients  o'  the  interval  ouadratic  equation  needed 
by  ROOTS. 

ROC'S  solves  the  interval  quadratic  equation  and  deletes  points  of  the 
Current  box  where  f(x)  >  F5AR. 

NEWTON  does  a  step  of  an  interval  Newton  metnod  to  find  stationary 
points  of  the  gradient  of  f (x )  in  the  current  box. 

I  IN’.  IF  is  an  IVSL  subroutine  used  to  compute  tne  approximate  inverse 
o'  the  center  of  the  Jacobian. 

UPDATE  evaluates  fix)  at  the  center  o'  each  new  box  and  updates  FBAR 
(the  smallest  known  value  of  f',x)). 

DEBUS  prints  various  information  useful  ir  debugging  the  proara"  or  <0 
evaluating  its  performance. 

PRINT!  is  used  by  DEBUG  to  write  appropriate  information. 

SECT!  is  called  by  ROOTS  to  obtain  the  intersection  of  two  intervals. 

SECT2  is  called  by  ROOTS  to  obtain  the  intersection  of  an  interval  with 
the  complement  (exterior)  of  an  interval. 

Note  that  each  subroutine  contains  comments  describing  its  function  and 

use. 
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Introduction 


This  report  I*  designed  to  serve  as  a  supplement  to: 

F.  D.  Crery  and  T.  D.  Ladner,  "A  Simple  Method  of  Adding  a  New 
Datatype  to  Fortran,"  The  University  of  Wisconsin,  Mathematics 
Research  Center.  U.S.  Army,  Technical  Summary  Report  *1065, 

May,  1970. 

A  description  is  given  here  of  the  changes  Mde  in  the  Fortran  pre¬ 
compiler,  CLUDGE,  when  it  was  implemented  for  the  IBM  3&0/o~  Washington 
State  University.  Throughout  this  report,  number*  in  square  brackets  refe 
to  page  numbers  in  the  report  mentioned  above. 


The  default  prefix  for  type  other  in  CLUDGE  has  been  changed  from 
000  (oh-2e ro-oh)  to  LNT.  [.5.  >5.  "33- 

Storace  Reou  ?  rement  f  or  T  vt»c  Other 

The  dcfau’t  storage  requirement  for  a  type  other  varab'e  has  beer 
changed  from  one  memory  word  to  two.  I&3* 

Copy  Character 

The  default  copy  character  has  been  changed  from,  a  s 'ash  (/)  to 

a  p^us  (-0.  D1-  *57- 

U-L®1  b’sed  bv  CIUDGE 

The  data  set  reference  number  for  SCRATCH2  has  been  changed  fror 
21  to  U  £l2,  loD. 

Sumr.a r y  P  r  1  r. t e d  by  CLUDGE 

The  summary  printed  out  by  CLUDGE  at  the  end  of  the  instructions 
to  CLUDGE  no  longer  inc’udes  section  one  of  the  symbol  tab'e.  C’O* 

Type  Statement 

The  default  type  statement  recognized  by  CLUDGE  for  the  new  data 
type  has  been  changed  from  "JTHER"  to  ihTRVAt.  C^3’ 


Pott i ng  Comments  O n  the  Instruct ions  to  CLUDGE 

The  character  which  $top$  the  examination  of  an  instruction  to 
ClUDGE  has  been  changed  from  the  do’ lar  sign  ($)  to  the  "at"  sign  {£) . 

p3.  233- 

JM 09  Con t r o 1  Cards 

The  statements  in  ClUDGE  which  process  H09  contro*  cards  have 
been  removed.  £’2,  ’>93* 

Octa 1  Cons  tants 

The  statements  in  CLUDGE  which  a’’ow  for  octa'  constants  have 
been  removed.  [2'3* 

Do’ lar  Sion  C per a  tor 

Operator  *7  has  been  changed  from  the  do’ lar  sign  ($)  to  the 

"at"  Sign  («) .  [33.  **6,  *9.  50.  59.  60.  6l ,  62.  66.  7^. 

Tne  S y*nt>o'  Tab  'e 

The  physica'  structure  of  the  symbo’  table  had  to  be  changed  since 
350  Fortran  allows  at  most  U  characters  to  be  stored  in  one  word,  not  6. 
Thus,  array  TABLE  in  ccxomon  block  SYMBOL  was  changed  from  a  '000  x  3 
f u 1 1  word  array  to  a  7  x  1000  halfword  array.  The  various  items  are 
now  stored  in  the  symbol  table  as  follows: 


\ 


Sect  to*  One  Ent r  ies  £36}: 


i 

>- 


Symbol  (2  chars  /  element) 


^  Type  Other  "Equivalent"  (2  chars  /  element) 


_  J 


Sect  ion  Two  Entries  ^3°3: 

Dimensioned  Variable  £3»3: 


Symbol  •< 

(2  chars  /  e  enent) 


V  - 


number  of  dimensions 


1st  dimension  (2  chars  !  element) 


Ccrnon  Block  Link  i 


<  cont inuat ion  link 


•  2nd  dimension 
(2  chars  /  element 


Nth  Continuation  Block  £373: 


B  lank 


Dimension  4n  -  2 
(2  chars  /  e lement) 


D  inens ion  kn  *  1 
(2  chars  /  e lenent)  * 


»  Dimension  4n  (2  chars  /  element) 


^  Dimension  4n  1  (2  chars  /  e'ement) 


(Continuation  Link 


Ar  1 1 tvne tic  Statement  Funct ion  £ 393: 


Symbol 

(2  chars  /  e lement) 
Type 

Number  of  Arguments 


A 

i 

I 

r 

4 

{ 


*  Expansion  Indicator 

•  End  Pointer 


Function  Subprogram 

r 

Symbol  \ 

(2  chars  /  element) 

( 

Type  { 


Result  Symbol 
(2  chars  /  e  .ement) 


Section  Three  and  Four  Entries 


M: 


j 


Beginning  Character  Position 


►  Type 

J 

1  Last  Character  Position 
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feature.  Other  routines  uslrg  this  feature,  DMC0DE  and  D“STCD  were 
deleted.  fll7j. 

i 1 08  FLD  Function 

To  facilitate  character  and  string  handling,  three  Assembler  routines, 
PUTC,  GETC,  and  KCmPAR  are  used  to  replace  the  1103  FLD  ■'’unction.  £ll8j. 


D  i  sc  la  tmer 

The  pert*  of  CIUDGE  pot  mentioned  ebove  »ere  left  etient lei ly  the 
.*  In  the  1I0S  vers  Ion.  Since  no  cl. In  -.*  9«-.n  to  the  correc.ne*. 
,he  I.OS  version,  no  gu.r.nt.e  cen  be  9iven  to  the  correctness  of 

the  360  version,  either.  Ho- . the  ch,n9.s  is.de  h.ve  been  thorov9hi¥ 

tested,  and  seem  to  be  working  correctly. 
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Introduction 


Anyone  who  has  ever  tried  to  do  numerical  calculations  on  a  computer 
is  probably  painfully  aware  that  his  results  are  at  best  only  appro* i-a- 
tions  to  real  solutions.  Many  types  of  error  generally  contribute  to 
this  phenomenon,  including  error  in  representet ion  of  Cota  within  the 
computer  and  error  inherent  in  the  way  the  computer  does  its  arithmetic, 
i.c.,  ‘'roundoff"  error.  Through  the  use  of  interval  ar  i  thmet  ic .  the  effe 
of  these  tv;o  types  of  error  nay  be  softened  somewhat.  Famon  Moore  [  ll 
give''  sene  interesting  details  on  the  use  of  interval  analysis  in  numer¬ 
ical  methods. 

The  basic  idea  behind  the  Interval  Arithmetic  Peerage  is  that  opera¬ 
tions  can  be  performed  on  interva  U  of  real  numbers ,  rather  than  on  the 
real  numbers  themselves.  These  operations  ore  done  in  such  a  way  that 
the  resulting  interval  contains  the  exact  solution,  and  (hopefully)  is 
the  smallest  machine  representable  interval  containing  the.  exact  solution 
which  can  be  calculated  from  the  data  given. 

This  report  explains  the  basic  capabilities  of  this  Interval  Arith¬ 
metic  Package  and  the  use  of  a  Fortran  precompiler,  CLUOCE  (pronounced 
with  a  long  "u")  which  allows  Fortran  variables  to  by  typed  "IhTP.w’Al" 
and  used  just  as  standard  variables  arc  used  in  Fortran  expressions. 
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I .  The  Concept  of  Interva 1  A r  i  t hme  t i c 

A r i thmet  i  c  oocrat ions  defined  on  intervals  are  merely  operations 
defined  on  the  two  sets  of  real  numbers  involved.  For  example,  addition 
of  two  intervals,  £a,b^  and  £c,d^,  is  the  set  of  sums  which  nay  be  ob¬ 
tained  by  adding  real  numbers  ,  x  and  y ,  where  a  £  x  b  and  c  «£  y  ^  d. 

That  is  to  say, 

[a.bj  ■»  [c.dj  “  <x  •»  y:  x  £  [a,bj  and  y  &  [c,d]>. 

In  a  similar  manner,  funct ions  nay  be  defined  on  intervals.  Thus 
for  a  real  valued  function  f,  f([a,b3)  is  the  interval  £c,dj  consisting 
of  all  points  y  such  that  y  ■  f  (x)  for  sane  x,  a  4-  x  f  b. 

The  set  notion  is  also  used  to  define  ccrpar isons  cf  intervals. 

For  example,  "e.bj  t  ^c,d^  means  that 

x  4.  y  for  all  x  £  [a,b]  for  a'l  Y  £  •e.d]. 

That  is,  must  be  true  for  all  possible  pairs,  x  and  y,  where  x  is 
taken  frcn  [*  ,bj  and  y  is  taken  from  Tc.dj. 

It  should  be  noted  that  an  interval  of  the  form  [a  ,a~  is  simply 
another  way  of  representing  the  real  number  a.  Such  an  interval  is  called 
deecnerate.  So,  in  a  sense,  interval  arithmetic  can  be  thought  of  as  an 
extension  of  real  arithmetic. 

I  I .  Add inp  a  New  Data  Type  to  Fortran 

Standard  Fortran  compilers  allow  programmers  to  use  several  different 
data  types,  such  as  REAL,  INTEGER,  COMPLEX,  and  DOUBLE  PRECISION,  when 
performing  numerical  ca  I  cu 'at  ions .  In  order  to  odd  a  new  data  type  to 
Fortran,  the  programmer  must  sonehow  simulate  operations  on  this  -ew 
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data  type  in  terms  of  existing  operations  and  data  types. 

A  convenient  way  of  adding  a  new  data  type  to  Fortran  is  to  use  a 
"precompi ler"  which  will  recognize  a  type  statement  for  the  new  data 
type,  and  then  translate  all  arithmetic  operations,  function  calls, 
ccxnpar isons ,  and  data  type  conversions  involving  the  new  data  type  into 
calls  to  the  appropriate  simulation  routines.  For  example,  a  precompiler 
would  accept  the  statements 

INTRVAL  A ,B ,C 


and  replace  them  with 


C  -  A  ♦  B 

COMPLEX  A.B.C 


CALL  InTADD  (A.B.C) 

where  INTADD  is  the  routine  which  performs  addition  on  the  intervals 
A  and  E,  and  stores  the  resulting  interval  in  C.  (T he  reason  for  re¬ 
placing  the  INTRVAL  type  statement  with  a  CPhPLEX  type  statement  will 
be  explained  later.)  The  Fortran  program  produced  by  the  precompiler 
Is  then  passed  to  the  Fortran  compiler  to  be  compiled  and  executed. 

The  routines  which  simulate  operations  on  Intervals  are  contained 
In  the  Interval  Arithmetic  Package,  and  are  explained  further  in  the 
following  sections.  A  precompiler,  named  CLUDGE,  is  available  which 
will  recognize  the  INTRVAL  type  statement  and  insert  calls  to  these 
Interval  arithmetic  routines.  Control  cards  and  the  deck  setup  for 
using  CLUDGE  are  given  in  Appendix  C. 

Note  that  routines  in  the  Interval  Arithmetic  Package  may  be  called 
explicitly  from  the  user's  program.  In  the  last  example,  the  user  could 


t  the  statements 


pu 

COMPLEX  A.B.C 
CALL  intaod  (a.b.c) 

directly  into  his  program,  and  then  a  precompiler  is  not  needed.  The 
precompiler  CLUOCE  has  been  implemented  only  to  make  the  use  of  the 
Interval  Arithmetic  Package  easier. 

1 .  The  Implenentat ion  of  Interva )  Ar  i  thnet ic 

The  following  sections  describe  the  methods  used  in  the  Interval 
Arithmetic  Package  to  implement  the  various  interval  operations.  The 
implementation  preserves  the  concept  of  interval  arithmetic  as  it  Is 
described  in  Sect  ion  I. 

Within  the  computer,  an  Interval  of  real  numbers  will  be  represented 
by  its  endpoints .  We  will  thus  refer  to  an  interval  [a,b3  as  an  ordered 
pair  of  machine  rcprcsentaole  single  precision  numbers,  a  and  b,  such 
that  a  £  b.  Hence  in  order  to  represent  an  interval  by  a  single  variable 
in  Fortran,  the  variable  name  rust  refer  to  two  memory  cells,  one  for  the 
left  endpoint  and  one  for  the  right  endpoint.  This  can  be  accomplished 
In  cither  of  two  ways:  by  declaring  the  variable  CCmPLEX  or  by  declaring 

*  S' 

It  to  be  a  real  array  of  two  elements,  i.e.,  COMPLEX  A  or  DIMENSION  A  ( 2) . 

If  the  user  is  calling  routines  in  the  Interval  Arithmetic  Package  directly, 
he  may  use  either  method.  If  he  is  using  the  preccrpi ler ,  CLUDGE,  he  must 
use  the  INTAVAL  type  statement  as  described  in  the  last  section. 

We  can  now  go  on  to  describe  the  various  interval  operations  In 
terms  of  the  representations  of  the  intervals  by  endpoints. 
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A  .  6 as  i_c  interval  A r  i  t  hme  t  i c 

The  four  basic  binary  operations,  addition,  subtraction,  multiplication 
and  division  are  defined  for  intervals  as  foltows: 

[a,b]  ♦  [c.d]  «[a  ♦  c.  b  ♦  d] 
ra,b  -  ^  c  ,d  *  a-d,  b-c"1 

[a,bj  x  [c.d]]  "[min  {ac  ,ad  ,bc  ,bd},  max  fac  ,ad  ,bc  ,bdj] 

[a,b[  /  Lc.dl  "Tmlnfa/c,  a/d,  b/c,  b/d_*  ,  rax  i  a/c  ,a/d  ,b/c  ,b/d; 

(note:  division  is  defined  only  if  0^£c,dj.) 

Note  that  these  definitions  are  consistent  with  the  genera'  set  definitions 

5 i ven  i n  Sect  ion  f , 

The  subroutines  in  the  Interval  Arithmetic  Package  which  perform 
these  operations  are  IN7ADD,  If.TSUB,  INTMUL ,  and  INTO  IV,  respectively. 
Appendix  A  describes  how  these  routines  may  be  called  directly.  If 
CLUDCE  is  being  used,  calls  to  these  routines  will  be  generated  whenever 
the  sy-bols  * ,  - ,  *  ,  or  /  are  encountered  and  cither  or  both  operands  are 
of  type  INTP.VAl.  (Note  that  mixed  node  _[_s  at  loved  oy  ClUDGE.  Convers  ions 
between  data  types  are  described  in  Section  111.  C.)  ClUDCE  ignores  a 
unary  *  sign,  and  ur.ar y  negation  is  performed  by  the  routine  INTNEG. 

For  further  information  on  the  theory  and  algorithms  used  for  this 
interval  arithmetic,  see  [lj  and  [2]. 

B  .  Bcs  t  boss  ibljs  F 1  oa  t  inc  Point  A r  i  t hmg  t  i c 

While  the  interval  operations  defined  in  the  last  section  are 
theoretically  correct,  they  cannot  simply  be  programed  using  the  computer’s 
binary  operations,  since  this  would  again  introduce  the  error  we  are  trying 
to  avoid.  Consider  the  following  exa~ple: 
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1.  _Rea_l  to  Inter  v  t  (INTC5B): 

A  degenerate  interval  is  created  whose  endpoints 
are  the  given  real  number. 

2.  Double  Precision  to  Interva  I  (INTC68): 

If  the  double  precision  number  can  be  expressed 
exactly  in  single  precision,  the  convers ion  is 
the  some  as  in  I,  If  not,  a  nondecencrate  interval 
is  obtained  by  rounding  dtxvn  on  the  double  precision 
number  for  the  left  endpoint,  and  rounding  up  for 
the  right  endpoint. 

3.  Integer  to  Interva  I  ( I h« T C A* 8 )  : 

The  integer  is  first  converted  to  double  precision, 
and  then,  to  an  interval  as  described  in  2. 

U ,  Corep  lex  to  Interva I  (INTC78): 

The  real  part  of  the  complex  nu-ber  is  converted 
to  an  interval  as  described  in  I.  The  imaginary 
part  is  ignored. 

5.  Interva  1  to  RcaJ^,  Doub  le  Prec  is  ic?"  .  or  I  ret  ece-  (lf.TC85, 

IHTC86,  ItlTCSfc); 

The  midpoint  of  the  interval  is  computed  in  double 
precision  and  is  then  converted  to  the  desired  type, 
rounded  if  necessary  to  single  precision. 

6.  Interva  1  to  Complex  (INTC87): 

The  real  part  is  set  to  the  rourded  value  of  the 
midpoint  of  the  interval,  and  the  imaginary  part 
is  set  to  zero. 


When  using  CLUDGE,  the  programmer  may  mix  variables  of  type  INTRVAL 
with  those  of  type  REAL ,  DPJBlE  PRECISION,  INTEGER,  or  CC-.PlEX.  Type 
IMTRVAL  dominates  the  other  data  types.  Thus  mixed  mode  arithmetic 
involving  a  variable  of  type  INTRVAl  will  be  done  in  interval  arithmetic. 
ClUDGE  will  generate  calls  to  the  various  conversion  routines  whenever 
necessary.  For  example,  suppose  R  is  real,  and  A  is  type  INTRVAL,  Then 
during  precompi lat ion  with  CLUDGE,  the  statement 

R  ■  2  *  A 


'rmrr 


is  replaced  by 
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REAL  1NTC8S 
C0MPLEX  A,  INTTEN (2) 

• 

• 

CALL  INTCA® (2 , INTTEM ( 1 ) ) 

CALL  INTKUL  ( INTTEM ( I ) ,A , INTTEM  (2) ) 

R  -  INTC85 ( INTTEM (2) ) 

Note  that  CLUDGE  will  allocate  temporary  storage  locations  in  the  array 
INTTEM  to  hold  intermediate  interval  results  whenever  needed. 

The  arithmetic  IF  statement  is  the  one  except  ion  to  the  fact  that 
the  INTRVAL  data  type  always  dominates.  If  the  expression  in  an  arith¬ 
metic  IF  statement  contains  a  variable  of  type  INTRVAL,  CLUDGE  will 
generate  the  statements  necessary  to  evaluate  the  expression,  just  as 
explained  above,  but  the  final  result  is  converted  to  type  REAL  and  the 
proper  branch  is  determined  on  the  br.sis  of  this  REAL  value.  For  example, 
for  the  statements 

INTRVAL  A 

IF  (A)  3.^.5 

CLUDGE  will  generate 

REAL  INTC85 
C0*‘.PLEX  A 
• 

IF  ( INTC85  (A) )  3,^,5 

The  user  ray  also  call  these  convers ion  routines  directly  from  his 
program  if  desired.  (See  Appendix  A). 

D.  Exponentiation 

There  are  four  routines  in  the  Interval  Arithmetic  Package  for 
raising  an  interval  to  a  power: 
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1.  Intcrva  1  to  Integer  Power  (INTXP**): 

The  result  is  obtained  by  successive  multiplications 
of  the  interval  with  itself.  However,  if  the  power 
is  an  even  integer,  the  left  endpoint  of  the  result 
will  never  be  less  than  zero.  For  example,  if  A  is 
the  Interval  [-1,2^,  then  the  result  of  squaring  A 
will  be  ^0,0,  even  though  A  multiplied  by  itself 
yields  £*2,AJ. 

2.  Interval  to  Rea  I .  Doub  le  Prec  is  ion  .  or  Ir.terva  I  Power 
(INTXP5,  IUTXP6,  INTXPv): 

The  result  of  A"  is  calculated  by  exp (Ex  log (A) ) 
where  the  interval  exponential  and  log  routines 
(described  in  Section  III.  F.)  are  used,  and  8 
is  converted  to  an  interval  beforehand,  if  necessary. 


3.  ReaJ,  Double  Prec is  ion ,  or  Inteoer  to  Interval  Power 
( INTXPS) : 


The  base  is  first  converted  to  an  interval,  and 
then  the  rout ine  INTXPS  is  used. 


The  precompiler  ClUDGE  will  process  occurrences  of  **  in  a  Fortran 
program  by  generat  ing  calls  to  the  appropriate  rout ines  ,  depending 
on  the  type  of  base  and  exponent.  If  CLUDGE  is  not  used,  the  user 
must  see  that  proper  conversions  are  done. 


E .  Cpmpar ison  of  Intervals 

The  Interval  Arithmetic  Package  contains  logical  functions  INTEQ, 
INTNE,  INTGT ,  INTGE,  INTLT,  INTLE  which  perform  the  comparisons 
.HE.,  .GT.,  .GE.,  .LT.,  and  .LE.  on  two  intervals.  These  relationships 
arc  defined  as  follows: 

£a,bj  m[c  ,dj  if  and  only  If  a*b»c*d 

fa,b]  t  rc.dl  and  °hly  If  b  4  c  or  «  s  d 

[a,b3  <Tc,d3  if  and  only  If  b  <.  C 

[a.bj  *  [c,d]  if  and  only  If  b  *  c 
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£a,bj  a.  [c,d}  M  and  only  If  a  *  d 

£a,b^  >  f  c.dj  if  and  only  if  a  a  d 

These  are  consistent  with  the  general  definitions  in  Section  I.  In 
particular,  note  that  if  R  is  any  one  of  the  above  relational  operators, 
then  [a,b7  R  £c,dj  dn{*  or,’y  'f  *  R  y  for  all  x  €  [a,b^  and  ye  fc.dj. 

For  example  ,  ifxtri.S^andyfi.^.1*}  it  is  not  dear  ho*  x  and  y  are 
related.  Hence  the  relational  operators  can  be  defined  only  on  disjoint 
interva's,  or  intervals  with  at  most  one  point  in  common.  Care  must  be 
taken  when  using  these  operators  since,  for  example,  ra.bj  <£c,d}  does 
not  imply  that  '"a.b]]  i  [\:,d}.  Even  more  confusing  is  an  expression  such  as 

. N?T.  (A.EQ.  E)  .AND.  .NCT.  (A.hE.B) 

which  is  t  rue  if  A  and  E  are  nondegenerate  intersecting  intervals,  even 
if  they  are  equal  a$  sets.' 

These  routines  may  be  referenced  directly  by  the  user  {see  Appendix  A), 
or,  they  will  be  referenced  by  CLUDGE  whenever  one  of  the  relationa’ 
operators  .EQ.,  .NE.,  .GT. ,  .GE.,  . LT. ,  or  ,L£.  is  encountered  with  a 
type  INTRVAl  operand.  {if  the  other  operand  is  not  of  type  INTRVAl,  it 
is  converted  before  the  comparison  is  made.) 

f.  Standard  functions 

The  Interval  Arithmetic  Package  provides  interval  counterparts  to  the 
Standard  Fortran  functions  SQRT,  EXP,  AS.PC ,  AlCGlO,  ARSIN,  ARCCS,  ATAN, 
SINH,  CCSH,  TANH,  CCS,  SIN,  TAN,  and  ATAN2.  These  routines  are  listed  in 
Appendix  A.  Thus,  if  desired,  the  user  nay  call  them  directly.  If  CLUDGE 
•  s  bemg  used,  calls  to  the  interval  counterparts  will  be  made  whenever  one 
of  the  above  routines  is  encountered  with  a  type  INTRVAL  argument. 


I 


s  12 

In  addition,  there  are  six  interval  routines  which  have  no  Fortran 

counterparts,  and  thus  will  never  be  referenced  by  CLUDGE.  They  are: 

INTSUP  Returns  the  right  endpoint  of  an  interval 

INTINF  Returns  the  left  endpoint  of  an  interval 

INTlEN  Returns  the  length  of  an  interval 

INTUPN  Returns  the  union  of  two  intervals: 

[a  ,bj  v  [c  ,d3  =  •  n  ( a  ,c)  ,  r.ax  (b  ,d)J 

INTSCT  Returns  the  intersection  of  two  intervals: 

[a.b^  n[c.d2  r  [ra«  (a  .c)  ,  min  (b,d'3 

(Tne  intersection  of  disjoint  intervals  is  not  defined) 

INTBND  Returns  an  interval  created  by  bounding  a  double 

precision  number  to  a  given  accuracy 

These  routines  are  also  listed  in  Appendix  A,  and  ray  be  called  directly 
by  the  user, 

(Note:  There  is  no  routine  in  the  Interval  Arithmetic  Package  which 
will  create  an  interval  from  two  real  numbers.  However,  this  may  be 
accomplished  as  follows: 

1.  Convert  each  real  number  to  a  degenerate  interval. 

2.  Form  the  union  of  these  two  intervals. 

To  create  an  interval  from  a  double  precision  number,  use  the  subroutine 
INTEND  and  specify  how  many  bits  of  the  double  precision  number  are  assumed 
to  be  accurate.) 

H.  Error  Conditions 

During  every  interval  operation  performed  by  a  routine  in  the  Interval 
Ar ,  t hret  ic  Package,  possible  error  conditions  are  monitored  and  a  flag  is 
set  whenever  an  error  occurs.  At  the  end  of  each  operation,  a  call  is 
r**Ct  to  INTRAP,  a  routine  in  the  Interval  Arithmetic  Package  which  takes 
•cpropr.ate  action  if  an  error  has  occurred.  Belcw  is  a  list  of  the  possible 


error  conditions,  alone  with  the  value  of  the  result  in  each  case,  and  the 
action  taken  by  INTRAP.  The  meaning  of  overflow  and  underflow  should  be 
clear.  An  "infinity"  error  occurs  when  an  attempt  is  made  to  round  up 
(down)  a  number  which  is  already  at  least  as  large  (small)  as  the  largest 
(smallest)  single  precision  number  possible.  The  symbo'  K  is  used  to 
represent  the  largest  possible  single  precision  number,  L  is  used  to 
represent  the  largest  possible  single  precision  number  which  can  be 
correctly  converted  to  an  integer,  and  C  indicates  a  correct  va'ue.  The 
possible  actions  taken  by  INTRAP  are: 

A  No  Act  ion  Taken 

D  Prints  Error  Message 

C  Prints  Error  Message  and  Traceback, 

and  I nc r e* en t s  the  Error  Courier 

Tne  fault  conditions  recognized  by  INTRA*>  are: 


Pounds  faults: 
Left  Endooint 


R  ioht  Endpoint 


Result 


Action 


no  fault 

infinity 

no  fault 

unoe  rf I ow 

overf low 

infinity 

inf i n  i  t  y 

no  fault 

infinity 

Overf low 

inf  in  i  ty 

inf  i n  i  t y 

infinity 

underf low 

under  f low 

no  fault 

underf low 

infinity 

underf low 

unde  rf low 

2.  Other  Faults 


Fault 


Result  Act  ion 


Division  by  Zero 
Zero  to  the  Zero  Power 
Square  Root  of  Negative  Number 
Log  of  Non-Positive  Number 
Underflow  during  Interval  to 
Rea  I  Convers ion 
Overflow  during  Interval  to 
Integer  Conversion 
Intersection  of  Disjoint 
Intervals 

ARCCS  or  ARSIN  Argument 
Out  of  Range 


undef ined 

[-K.K] 

[-K.»0 


L 

undef i ned 


C-k.k] 


c 

B 

C 

c 

A 

C 

c 

c 


Provisions  have  also  been  made  to  allow  the  user  some  control  over  the 
action  taken  when  errors  occur.  These  are  explained  in  more  detail  in 
Appendix  D. 


I.  tnnut  and  Output 


The  Interval  Arithmetic  PaCKaoe  contains  no  special  routines  for 
inpjt  and  output  of  interval  variables,  so  the  user  must  devise  his  own 
methods.  This  is  a  stra  ightforward  task  and  no  more  difficult  than  the 
task  of  input  and  output  for  ordinary  Fortran  variables  as  long  as  the 
user  remembers  two  things: 

1.  Interval  variables  require  two  memory  cells  and  hence 
must  be  treated  as  either  complex  variables  or  arrays  of 
two  elements.  If  they  are  declared  as  arrays,  they  may 

be  read  in  or  written  out  ju-St  as  ordinary  arrays  are.  If 
they  are  declared  COMPLEX,  they  may  be  read  in  or  written 
out  just  as  complex  variables  are.  If  CLUDGE  is  used,  all 
variables  of  tvne  INTRVAL  are  automatically  declared  COMPLEX 
and  should  be  treated  so  in  all  input  and  output  statements. 

2.  Fortran  FPAMAT  statements  allow  the  user  to  specify  the  number 
of  digits  to  be  written  out,  and  so  what  is  actually  written 

•  s-  the  value  rojnded  to  that  many  places .  Care  should  be 
taken  to  insure  that  this  feature  does  not  unintent ional ly 
misrepresent  the  actual  endpoints  of  the  interval  computed. 
That  is,  rounding  up  on  the  left  endpoint  or  down  on  the 


right  endpoint  could  cause  the  interval  to  appear  smaller 
than  i t  actua 1  ly  is . 

Note: 

Arty  appearance  of  a  variable  of  type  INTRVAL  in  a  READ  or  WRITE 
statement  causes  the  message 

*****  WARNING  -  TYPE  OTHER  ARGUMENT  IN  I/O  LIST  ***** 
to  be  printed.  This  message  is  simp'y  a  reminder  to  the  user  that  a 
type  INTRVAL  variable  is  being  used  and  must  have  two  fields  reserved 
for  it  in  the  FORMAT  statement.  The  message  ooes  not  indicate  an  error 
in  syntax,  but  merely  points  out  a  frequent  source  of  l/C  errors. 
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from  the  Washington 
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rwe»{«  Equivalent  if  feet 


Store  INTSTR  Subroutine  (AI.ARES)  '  •  *ARES  -  Al 

I  t 

Negate  '  INTIlEG  Subroutine  (AI.ARES)  !  -  (unary)  JARES  •  -Al 
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A ppend ix  B  * 


Jnterva!  Library  Routines 


In  the  following  table,  the  values  in  the  range 
approximate  extreme  values  of  the  interval  function, 
raximum  single  precision  floating  point  number.  The 


column  are  the 
M  stands  for  the 
number  in  the  "Errors" 


column  is  the  value  of  the  fault  flag  for  this  error  condition. 


Name 

1 

Range 

Error 

Explanation  or  Other  Comments 

Library 
Rout i nes 

‘Assumed  Ac 
racy  in  #  « 

IhTACS 

[O.tt] 

23 

Argument  out  of  range:  -1  to  1 

Used 

DARCPS 

b  i  ts 

52 

intass 

23 

Arcu-ert  out  of  range:  -1  to  1 

OARS  It. 

52 

|!.TA  Tli 

-- 

— 

DATAN 

52 

If.  TAT  2 

v 

16 

Division  by  jero  attempted 

DATA!! 

52 

INTCC?  i 

t- 1 . 0  ■ 

-- 

— 

DEC’S 

AS 

INTCSM 

f  i.«‘ 

-- 

1 f  aas.  va l .  of  endpoint ^  175.0, 

DCPSH 

51 

INTEXP 

[o.rt] 

175.0  used 

If  abs.  val.  of  endpoint  *  170.0, 

DEXP 

52 

INlLN 

[-181. 

'9 

170.0  used 

Log  of  non-positive  nu-  er 

DlOG 

50 

INTlPS 

’  175] 

[-79.76] 

19 

Log  of  non-positive  number 

dlcg  1 0 

except  25  < 
f.  99?  5i  1.00C 
50 

ll.TSIN 

C->.0 

DECS 

except  25  • 

r.9?95.  1.00C 

AS 

If.TSNH 

l-M.H] 

-- 

If  abs.  val.  of  endpo  int  »  1 75. 0, 

OS  IHH 

51 

INTSQT 

:  o.m] 

19 

175. 0  used 

Square  root  of  a  negative  number 

DSQRT 

except  A0  i 

f-.ooos,  .0: 
53 

If.TTAf. 

10 

Error  when  interval  argument  spans 

DTAt; 

33 

INTTfiH 

1 

c-«.o 

the  discontinuity  of  tan  at  -?$  or 
at  v/3  .  The  standard  function 

DTAN  will  also  generate  an  error 
if  either  endooint  is  "too  close" 
to or  to  v/i  .  No  check  for  this 
is  made  in  INTTAN. 

0TANH 

52 

Hcxv  to  Use  CLUDGE 


The  Fortran  prccomp  i  !cr ,  CLUDGE,  will  accept  almost  aM  Fortran 
i  tatr-cr.ts  allowed  by  the  IBM  Fortran-G  compiler.  The  exceptions  are 
c iven  here: 

1.  CLUDGE  will  not  accept  any  of  the  following  statements: 

JMTEGER-2 
INTEGER-  ** 

REAL-  *4 

real-3 
Cti  t?LEX  S 
CP'.PL  ’  6 
ICG  1  CAL-- 1 
LOGICAL  A 

CLUDGE  does,  howt-ver,  accept  the  following  type  statements: 

REAL 

INTEGER 

CP-mPlEX 

BP  JalE  PRECISION 

LOGICAL 

IN7RVAL 

Thus  the  user  shou’d  avoid  us  inc  INTEGER  2,  L0G  ICAL-'  I , 
and  COMPLEX*! 6  data  types  a’together  when  using  CLUDGE,  and 
should  replace  other  '*■.•*'  type  statements  by  ordinary  type 
statements,  i.e.,  use  REAL  rathrr  than  REAL'-A,  etc. 

2.  The  following  statements  arc  accepted  by  CLUDGE,  but  no  further 
act i oh  is  toLcn  on  them; 

DATA 

equivalence 

external 

FORMAT 

namelist 
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In  particular,  extreme  caution  should  be  taken  if  type  IKTfU'AL 
variables  are  used  in  DATA  or  EQ.L*  I  VALENCE  statements,  since  the 
form  of  these  variables  is  altered  by  CLUDGE.  Also,  since  F PR", AT 
statements  ore  recognized  but  not  processed  further,  the  user 
Shou'd  ovoid  using  FDR’tAT  os  r  he  name  of  on  array. 

3.  The  main  program  must  be  the  first  input  to  CIUDGE,  and  at  most 
one  main  p-oerc-  roy  appear  in  the  deck. 

li.  The  ra-ics  INTTE“  and  I!  7P.LT  arc  generated  by  CIUDGE  for  special 
purposes,  and  should  not  be  explicitly  used  by  the  progra-Tner . 

5.  If  CL.DGE  must  generate  additional  statement  numbers,  they  will 
start  with  or  the  first  integer  a'tcr  30000  which  has  not 

already  been  used.  Thus  the  programmer  should  attempt  to  use 
statement  numbers  belo.v  30000  whenever  fcssiblc. 

CtUDGE  does  a  minimal  amount  of  error  checking,  but  may  occasionally 
produce  error  messages.  These  messages  ere  listed  in  Appendix  E.  Note 
that  the  type  OTHER  referred  to  by  these  messages  is  type  INTRVAL. 

ClUDGE  is  actually  a  neneral  purpose  precompiler  for  FORTRAN,  and 
accents  several  instructions  which  will  alter  the  way  in  which  it  procestes 
r  program.  These  instructions  arc  described  in  [$j  and  [ 7 j .  but  the  inexperi- 
in.-td  programmer  should  not  attempt  to  use  them. 


setup  for  CLUDGE  is  given  below 


0™  ft  if 
nr  data 
cares 


-  Data  Cards 


Source  Cards 


II  INTERVAL. SY$ INaDDi- 


//S  tcpr.o-ic.:E>  E  C_  INTERVAL 


x ; pb  card 


.  .  I  #»  a  r  A 


Deck  Setup  to  Prcccnpile,  Compile,  and  G? 


S  tcrr..vnc 


Stuff c c  Cards 


.a^V 


?k'jt 


may  be  any  name  of  I  to  8  characters  (first  character 
must  be  alphabetic)  or  r.a y  be  left  blank. 


denotes  at  least  one  blank. 


arc  the  Fortran  program  which  uses  the  IMTRVAL  type 
s  tat cneni . 


A  program  run  under  this  control  card  setup  will  produce  three  levels 


ef  output: 


| .  Put  pi- 1  f  rrT  C  ll'PCS  : 

This  is  o  listing  of  the  source  deck,  along  with  any  error 
messages  which  ClUCGE  generates. 

2.  Output  f  rc-  the  Fortran  C  p<~-n  i  lr; 

This  is  a  listing  c‘  t*-c  Fortran  program  which  CLUDGE  produces 
by  changing  the  type  luTRVAl  statement  to  CC'.FlEX  and  replacing 
opera:  ic~s  on  type  INTF.VAi,  variaples  with  subroutine  ca'ls  s-'d 
function  rc^cre'ces.  Any  erre'  messages  produced  by  me  Fortran 
compiler  will  be  riven  here. 

3*  Cutout  f  rev  E>e;ut  I  on  t 1  c  Prorrar; 

This  is  the  cutout  which  is  generated  by  the  execution  of  the 
compiled  Fortran  prcg-a~.  Any  error  messages  generated  by  this 
execution  will  be  given  here. 


These  three  levels  of  output  are  illustrated  in  the  sa~pie  pregron 
which  follows. 
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Source  iiiii'f  k/if(  ItHff  0f  orioi  wte i Safes  +JJcJ 
^  CLI/OC  t. 


1  B*'  3fc0/fe7  version  oc  cludde 


S£W plf  PROG0*'*  to  C^MPJTF  A  TABLE  OF  FACTORIALS  AND  THEIR  MAT  UR 
LOGARITHMS  US1MS  INTERNAL  ARITHMETIC 

:f  al  1ST  P-'F,  IMTSUP.L'FT 
|  1  3  V  a  L  FACT, t*:F  ACT 

%.*  \  Tf  (  ,  1  > 

‘ AC  T«  1  . 

D?  10  1*1.25 

rfci*FACT*i  Jif  p,  if  \ot  fyp1a#tti'on 

i  ‘if  act*  ai  or.  i  fact  )  r  • 

l  t  F  T*  I \ T  I  *  F ( c  ACT  J 
«  1  f.H’  *  1  NT  SLR<  F  ACT  I 
) r  | Tf ( f ,2  I  I , L c FT . R 1CH T .IMF  ACT 

•••••  wAPMING  -  TVPE  OTHFF  ARGUMENT  I\  1/0  LIST  ***t* 

\0TF  THAT  ThE  PRINTING  C>F  THE  VALUES  OF  THC  TUT  Tvor  JnT»VAL 
VARIABLES  IS  PTOf.-'  *Z  D  IM  OlFFEPfNT  *AVS 

1  r  lb-AT  I  •  I  T  AP.L  E  OF  INTERVAL  FACTORIALS  AMD  THEIR  L  0" &°1  THMS  •  ,  /  /  /  , 
I  '  M‘,16X,*M  f  ACTOR  I  AL  '  ,  22X  ,  'l  0G(  N  FACTOR |AL I  » ,// I 

:  (  oa»*at  (  '  •  .1? .?< <.x,m  •  .E is. s, • , •  .eis.e. •  i  •  n 

CALL  F  X  I T 

t  \D 


?"]***  f rciuetd  by  UUDZt  f„  ^ 

foATfiAV 


i«vft 


MAIN  DATE  =71166  15/2.' 


C 


( 

f 

C 


c 

r 

i 

r 

( 


•  ••••  PROCESSED  PY  1°m  060/67  VERSION  OF  Cl  UDGE  *•••• 

C  P  m  p  i  E  x  I  N I  T  E  v  (  2  ) 

SEAL  !  \  T  1  N  F  «  1NTSJP.  LEFT 
C^wPlf*  PACT,  LNFACT 

SAMPLE  o^rrpt^  T  0  COMPUTE  A  T£«UE  PF  FACTORIALS  A\D  THEJR  nATlT 

LOSA3IThvS  USING  INTERVAL  ARITHMETIC 

*o ITF( 6, 1 » 

CALL  I  NT  C  5  ?  (1.,  FACT) 

PP  10  1*1  ,?S 

CALL  1  \  T  C  a  3  (1,  INTTEMIIH 

CALI  INTMUL  (FACT,  !  NTT  EMU),  JNTTFMI  211 

CALL  I  \T  S  T  p  (If.TT£M(?),  FACT) 

FALL  1  Mln  (FACT,  L’.FACT  I 
Itri  *  1NTINF  (FACT) 

?  1  F.sT  *  1  NT  SUP  (FACT) 

K  V  P  I  T€(6,?)  1  ,LEFT,S  jr.HT  ,LNf  ACT 

•  ••«-*  WARNING  -  Tvof  CT^fr  ARGUE  NT  IN  I/O  LIST  ***** 

NOTE  THAT  THE  PRINTING  OF  THE  VALUES  0s  THE  TwC  TYPE  INTRVAL 
V  AR I  A  8L  c  S  IS  PROGRAMMED  IN  DIFFERENT  WAYS 


1  rfjMMAT  (  •  1  T  A«LE  of  INTERVAL  FACTORIALS  and  THEIR  LOGARITHMS 
•*  N  '  ,  1  6  x  ,  ’  N  FACTOR  I  AL  •  .22«»  •LO'il ‘J  T  AC  T  pc  1  ALT  »  ,  /  /  ) 

2  FO«imAt  (•  •  ,  !2,2(*x,  '  (  '  ,E  ,  •  »F1S.  9,  •)  '  )) 

CALL  -  *  I T 

END 


VtT 


///, 
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Outfit  f'oyam  execution. 


■;  intdVAL  FACTORIALS  AMD  Th=I»  LOGARITHMS 


N  factorial  LOGIN  factorial) 


r.icr  c^cf 

n  , 

C  .  1 0CD00O0E 

01  ) 

1 

C.  0 

9 

O.o  ' 

) 

r.vcic  cgcf 

'U  . 

0.2903C0C  E 

Cl  ) 

l 

C.6R7  1  A-Fl?c 

CD, 

D.6D31A71PF 

DC) 

F  .  *■  *  T  f  DDL  OF 

01 . 

o.sonoooo  1 

01  ) 

( 

C.  175  1  ’5?Fr 

01 , 

0. 170175C5F 

Cl  ) 

C  DOOOF 

0?  . 

D.  ?  <.0D'  TO0F 

:? » 

( 

C.  31  7=  '.5?«F 

Pit 

0.31  7  6  0  5  3  c  r 

C  1 ) 

**.  1  ?:::<■  d?c 

0  3  1 

0.1  2C DOOGDF 

03 ) 

I 

0. A7F  7  A  DDP  r 

:  1 , 

C.A7e74c;fr 

Cl) 

r  7  i  f  r»'  '^''C 

03  • 

C  .  72C TCrD''F 

031 

( 

0.  55  ’S  ?  5 "3  F 

Cl  , 

0.5579251  3  = 

Cl  ) 

r  ref 

OA , 

D.5  D-OODfCF 

DA  » 

1 

C. F5?F 1 60=  = 

01 , 

0.65251617= 

Cl) 

f..*f  ??<_C  CDF 

05, 

0.401230CCE 

05  ) 

( 

G. )D5'A6C3E 

02, 

D.  1C6DA«,0a  = 

C  2 ) 

f.'PHl  CDE 

Cs. 

'. 362S*::CF 

05  ) 

1 

0.  12FC 1 327C 

C  7  , 

0.  12SD1 52  Fr 

02) 

F-.  >f  ?  63'  CDF 

07  , 

0. 362B80GDF 

0  7  I 

1 

0. 15 1: aa) rr 

r,  •> 

•  *  9 

D.  1  SlCAAl’E 

02) 

r.  '  =  •♦  1 1  roo  = 

O'* , 

G. *9916900* 

C°  ) 

< 

0.  1  7 CD  2  3'' a  r 

C  ?• 

9.|T90?3U*f 

02) 

' . A?:  'C  1  60F 

. 

D • A TQDDIACE 

OR) 

I 

0.199172 | 3E 

•  n 

V  F  | 

0. 199=~22e= 

02) 

r . * ’? 7C 1 ?7£ 

ID. 

C.O??702 1 F= 

10) 

( 

0.7?rc2155= 

02, 

:. 2255217:= 

0  2 ) 

6  7  1  7R2  1  6f 

1 1 . 

D.P717P7A6E 

1  1  1 

( 

0.2515120°” 

0  7, 

2.25191223= 

02) 

C  .)  ?*  707?PE 

13  . 

:■  .  1  3C  7fe7  *>  f  c 

131 

< 

0.27FRR261F 

C  2 , 

D.  27porl277  = 

C  2  ) 

a.?;-'*??765E 

1*., 

0.2DR22P16F 

1  A  1 

( 

0.  3D 5 7 1 PA4C 

02, 

0. 30671 S75F 

02  1 

r.nt..,,  o*.qif 

)F>, 

0. 3555t>7'3<3t 

15) 

( 

C.  33505766  = 

02, 

2 . ?75D50F1= 

02) 

'  A*  ?3hl  7f 

It. 

C.OA023R7SE 

1  6  1 

( 

C.  35  3=  CA7  ’  = 

02, 

0. 3<63°5A62F 

02) 

r.  1  (• )  F 

1  F  . 

D.  1  21  6A536E 

16) 

( 

C.3R33c57aF 

02, 

D  .  3R33'aoqr>  = 

02) 

'  .  ?A3?^'»i>rE 

10, 

C.?A3?R0P1  = 

1  R  1 

( 

0. A2’3  56'?  = 

02, 

D.A2335632F 

0?) 

<  .  M  ;  a  D  h  f'a  = 

2-, 

0.5 1CR1GP6F 

2D) 

( 

O.A53*r 1?7' 

02, 

P.A53801<-2" 

02) 

*  .  II  ??'/<•  77f 

22, 

'■'. 2  1  240'3  =>= 

22  ) 

( 

0  •  A  5  A  7  l  1  76  = 

02, 

0.A5A71  1=>1  = 

f  2) 

•'  51  R45E 

23, 

C.25ES?ooor 

23) 

( 

0. 51 5  3  065°  c 

0  2, 

D .  51 60666°: 

0?) 

f  .*  2  •><.<.#  *,«,= 

2A  , 

0. 6 2 Da 501 «  = 

2A) 

1 

G. 5A76 A71 AF 

02, 

D. 5A7?A7tA= 

02) 

M’Mll  56E 

20, 

0. 1  3  51  1  2  5SE 

26) 

( 

0.56CD3601F 

C  2 , 

C  .5  =  00361  62 

C  2 ) 

1** 

g**  »***v 

V— 


"  - - - 


1 
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An->e*di*  P: 


Error  C . ntrol 


Eelow  is  a  complete  table  of  the  faults  which  are  recognized  by 
IKTRAP  The  value  of  the  result  after  each  fault  has  already  been 
•Sven  in  Section  III.H. 


|  Fault  flag 

1 

Fau  1 1 

.  _ 

Cond i t ion 

Default 

Response 

left  endpoint 

right  endpoint 

0 

no  fault 

no  fault 

m  * 

1 

no  fault 

overf low 

4* 

2 

no  fault 

infinity 

3 

3 

no  fault 

underf low 

0 

4 

overf low 

no  fault 

4* 

5 

overflow 

overf low 

4* 

6 

ove  r  f  1  cm 

infinity 

3 

7 

overf low 

underf 1 ow 

4*.  j 

8 

-infinity 

•  •  *  —  •  • 

no  fault 

3 

9 

-inf i n i t  y 

overf low 

3 

10 

-infinity 

infinity 

3 

1 1 

-  inf  in i ty 

underf I ow 

3 

12 

unde  rf 1 ow 

- 

no  fault 

0 

13 

unde  rf 1 ow 

overf low 

4* 

14 

unde  rf low 

infinity 

3 

15 

underf low 

underf low 

0 

•  l 6 

d i v is i on  by  zero 

3 

17 

zero  to  the  zero  power 

1 

18 

square  root  of  a  negative  number 

3 

19 

log  of  a  non-positive  number 

3 

20 

underflow  during  interva 1 -to- rea 1 

0 

1  21 

overflow  during  interva 1 -to- inteqer 

3 

22 

intersection  of  disioint  intervals 

3 

i  23 

ARCCS  or  ARSIN  argument  out  of  range 

3 

*  cannot  logically  occur 


0-2 


The  action  taken  by  INTRAP  is  determined  by  the  value  of  the  array 
H0NT0R  corresponding  to  the  value  of  the  fault  flag: 

MONTORdFAULT)  Actio_n 

0  Exit 

1  Print  error  message 

2  Print  error  message  and  traceback 

3  Print  error  message,  traceback, 

and  step  error  counter 

1«  Print  error  message,  traceback, 

and  stop 

The  user  may  communicate  directly  with  ll.TRA®  through  comm, on  block 

INTTLT.  The  statement  needed  is 

CCfrtHON/ INTFlT/IFAUlT, NAME, PARA'lS  (2,3) ,  |PAR,PAR,DPAR,MONTOR(23)  , 

«  KNTR  (23)  ,K0iJi«TS  (23) 

where  IFAULT  is  the  value  of  the  fault  flag,  NAME  indicates  the  routine 
in  which  the  error  occurred  (see  Appendix  A  for  the  va'ues  of  NAME),  the 
arrav  PARA'iS  contains  the  interval  arguments  (up  to  three  of  them,  stored 
es  columns),  I  PAR  contains  the  possible  integer  argument,  PAR  the  possible 
real  argument,  and  DPAR  the  possible  double  precision  argument.  Obviously, 
the  type  of  arguments  oassed  depends  on  the  type  of  arguments  of  the  rou¬ 
tine  in  which  the  error  occurred.  M0NT0R  gives  the  default  response  for 
each  of  the  possible  23  errors.  The  values  in  this  array  may  be  altered 
dynamically  by  the  user  to  change  the  action  taken  by  INTRAP.  For  example, 

M0NT0R(l6)  -  1 

causes  INTRAP  to  simply  print  an  error  message  when  division  by  zero 
(fault  "16)  is  attempted.  The  array  KNTR  gives  the  maximum  number  of 
tires  each  error  may  occur  before  execution  is  terminated.  This  array  may 
••so  be  altered  dynamically.  Thus  the  statement 

KNTR (IS)  •  U 

allots  taking  the  square  root  of  a  negative  number  to  be  attempted  four 
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i inei  before  the  iob  is  Hatted.  The  default  value  for  this  array  is  1 
for  each  error.  The  array  K0UNTS  contains  the  number  of  times  each 
error  has  occurred  during  execution.  Whenever  the  cor  respond ing  value 
of  MONTOR  '*  3  for  a  given  error,  and  that  error  occurs,  then  the  appro¬ 
priate  value  in  KPJ’tTS  is  incremented  and  tested  against  the  one  in  KNTR. 
if  the  maximum  numper  of  occurrences  for  that  error  has  been  reached,  the 
}ob  •'  terminated  by  It.TRAP. 

hence  if  the  default  opt  ions  for  INTRA?  are  used,  any  error  having 
a  corresponding  va'ue  of  3  or  U  in  MONTOR  will  be  fata’. 


Error  and  D  iaonos t  i c  Messages  from  CLUDGE 

This  append i *  gives  a  list  of  all  the  error  and  diagnostic  messages 
Included  ih  the  CLUDGE .  They  are  listed  alphabetically  according  to  the 
first  invariant  information  contained  in  the  message.  Variables  are  indi¬ 
cated  by  xxx  in  the  case  of  na-es  and  numbers  (occasionally  M  and  I.  are 
kied  for  numbers),  end  X  or  Y  in  the  case  of  characters.  Ue  a’so  give  the 
-eaning  of  the  message  if  it  is  not  se If -explanatory  and  the  meaning  of  any 
variab’e  parts  of  the  message,  further  explanations  cr  the  table  names  and 
code  words  used  in  the  messages  car.  be  found  in  [$1.  The  wores  TYF£  CTr.ER 
In  the  messages  refer  to  TYPE  1NTRVAL. 

ARITHMETIC  STATEMENT  FUNCTION  XXX  REFERENCED  VlTh  I.  ARGUMENTS  -- 
M  AROU-EhTS  DECLARED 

attempt  to  pop  empty  operator  stack 

Atte-pt  to  compile  a  Syntactically  improper  statement 

attempt  to  pop  empty  operand  stack 

Attempt  to  compile  a  Syntactically  improper  statement 

ATTEMPT  TO  READ  FAST  END  OF  FILE 
Pre-compi ler  error 

CC-PiiER  ERROR  ILLEGAL  OPERATOR  IN  STACK 

(1)  Closing  parenthesis  ')'  in  steel: 

(2)  Operator  associated  with  commas  is  not  a  function  call 

COMPILER  ERROR  <rw  LOOP  IN  COMMON  BLOCK  LINKS 

Common  block  links  have  been  improperly  constructed  resulting  in  a 
circle  of  entries  in  some  common  block 

COMPILER  ERROR  **  $  OR  )  IN  STANDARD  CODE  LIST 
Code  list  was  improperly  generated 

CONSTRUCTION  ERROR  --  UNEXPECTED  LETTER  FOLLOWS  CONSTANT 

The  character  following  a  string  interpreted  as  a  constant  is  not 
an  operator 
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CONVERSION  WITH  TYPE  LOGICAL  UNDEFINED 

One  side  of  an  arithmetic  statement  function  declaration  is  type 
logica1  and  the  other  has  a  standard  data  type  (other  than 
logical) 

DO  LOOP  NOT  ENDED  OR  BADLY  NESTED 

•X'  DOES  NOT  FOLLOW  ' Y '  IN  ALPHABET 

A  construct  of  t ne  form  (Y-X)  was  found  in  an  implicit  type 
declaration  where  'X'  pr.  eeds  ‘Y*  alphabetically 

END  INSERTED 

Warning  message  *  -  a  new  program  unit  was  begun,  but  no  end  line 

was  present  in  the  proceeding  unit  and  no  '•••NO  ENOS'  instruction 
was  present  in  the  instructions  to  the  CLUOGE 

ERROR  IN  CONSTRUCTION 

Generalized  error  message  indicating  error  in  the  construction  of 
a  constant 

ERROR  IN  INSTRUCTIONS  TO  PRE-COMP  I LER . 

EXECUTION  DELETED, 

Fata!  errors  were  encountered  in  the  instructions  to  the  CLUDGE 

HOLLERITH  CONSTANT  EXTENDS  PAST  END  OF  CARD  IMAGE 

A  Symbo’  interpreted  as  a  Hollerith  constant  has  too  great  a 
character  count 

HOLLERITH  CONSTANTS  MAY  NOT  AP»£AR  IN  STATEMENT  FUNCTIONS 

I  THINK  you  FORGOT  SO-ETHING 

No  data  has  bee-  pro. iced  for  the  CLUDGE 

ILLEGAL  COMBINATION  CF  OPERATORS 

An  il'egal  combination  of  operators  has  been  found  --  see  the 
description  of  the  parsing  method  and  decision  table 

ILLEGAL  COMMON  CLOCK  name 

Common  block  r.p-e  is  not  blank,  an  integer  constant,  or  an  identifier 

ILLEGAL  LIST  ITEM 

A  dimension  or  type  dec'arat'on  statement  is  improperly  constructed 
or  punctuated 

ILLEGAL  TYPE  NAME 

An  improper  type  name  has  been  specified  in  a  *type  instruction 
(note:  The  test  general ing  this  message  is  not  exhaustive) 

illegal  variable 

An  identifier  begins  with  a  character  other  then  a  letter 

improper  appearance  cf  decimal  point 

A  decimal  point  or  period  was  encountered  unexpectedly  whi'e  examining 
a  constant 


IMPROPER  APPEARANCE  OF  UNARY  OPERATOR 

An  operator  (other  than  '(')  immediately  follows  a  or 

IMPROPER  FILE  ASSIGNMENT 

Two  files,  other  than  output  and  print,  have  been  assigned  to  the 
same  logical  unit 

INVA'.IO  OPERATOR  (XX)  ENCOUNTERED  IN  TYPE  OTHER  CODE 

Pre-co-pi  ler  error  --  a  parenthesis,  comma,  logica',  or  relational 
operator  appears  with  a  type  OTHER  result  indicated.  XX  is 
the  operator  number  found 

LOGICAL  if  STATEMENT  FOLLOWS  LOGICAL  IF  STATEMENT 

XXX  MAY  NOT  BE  DELETED 

The  intrinsic  function  named  XXX  may  not  be  deleted  by  the  user 

MORE  THAN  ic  CONTINUATION  CARDS 

NEW  INTRINSIC  INCOMPLETELY  SPECIFIED 

NON-TERMINATING  DIMENSION  SIZES 

Tnc  list  of  dimension  sires  is  not  foi'owed  by  a  ') ’ 

XXX  NCT  IN  INTRINSIC  TABLE 

Function  XXX  specified  in  a  delete  instruction  cannot  be  found 
in  intrinsic  tab'e  or  is  a  type  OTHER  intrinsic  function 
(non-f  a ta 1 ) 

OVERFLOW  CF  CODE  TABLE  --  XXX  lines 

Not  enough  roar  in  the  code  tab'e.  XXX  is  the  va'ue  o*  CCDS  I Z 

CVERPlOJ  OF  CORRESPONDENCE  TABLE 

Not  enough  tab'e  space  to  associate  a  user  defined  statement  'abe< 
terminating  a  'DC1  ’oop  with  a  temporary  interna'  statement 
number  (non-fatal  unless  caused  by  a  statement  with  a  copy 
character) 

OVERFLOW  OF  DO  STACK  WITH  n  ENTRIES 

DO  loops  are  nested  too  deeply.  'N ’  is  the  depth  of  nesting  allowed 
by  the  CLUDGE 

OVERFLOW  IMAGE  --  XXX  CHARACTERS 

Extended  area  of  image  exhausted  while  adding  to  section  U  of  the 
symbol  table.  XXX  is  the  value  of  IMS 12 

OVERFLOW  OF  OPERATOR-OPERAND  STACK 

Space  for  the  combined  operator-operand  stack  has  been  exhausted 

OVERFLOW  OF  STATEMENT  FUNCTION  CODE  TABLE.  XXX  ENTRIES 

Skeleton  of  expended  arithmetic  statement  function  cannot  be  stored 
in  skeleton  table.  XXX  is  the  value  of  STSIZ 
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OVERFtCW  OF  STATEMENT  NUMBER  TABLE 

More  than  'RNGSIZ'  user  defined  statement  numbers  have  been  declared 
in  the  ranee  that  CLUDGE  must  remember  (r, on-fatal,  but  may  be 
the  cause  of  duplicate  statement  numbers) 

OVERFLOW  OF  SYMBOL  TABLE  --  XXX  ENTRIES 

Not  enough  symbol  tab’e  space  for  a  new  entry.  XXX  is  the  value 
of  TAB, ".AX 

PASS  2  END-CF-FILE  ON  SCRATCH  2 

Pre-conpiler  error  or  program  unit  with  no  correct  statements 

PASS  2  EQUIVALENCE  TA?lE  OVERFLOW 

Not  enough  table  space  to  store  equivalences  between  internal  and 
external  statement  ’atel  numbers 

PASS  2  S-FlAC  PRECEEDS  D-FLAG  ON  SCRATCH  1 
Pre-comp i ler  error 

SYMBOL  PREVIOUSLY  ASSIGNED  TO  A  COMMON  BLOCK 

A  synbo'  has  teen  declared  to  be  in  COMMON  more  than  once  in  the 
same  program  unit 

SYNTAX  ERROR  IN  ARGU"EnT  LIST 

Tr.e  argument  list  in  the  dec'aration  of  an  arithmetic  statement 
function  is  not  properly  punctuated 

TOC  MANY  ARGUMENTS  (XXX)  FOR  TABlE 

The  argument  list  in  an  arithmetic  statement  function  is  too  'one 
to  fit  in  the  tab'e.  XXX  is  t^-e  .maximum  number  of  arguments 
that  will  fit  irto  the  table. 

TOC  MANY  CONTINUATION  CAROS  --  *£n0  POSSIBLY  MISSING 

TYPE  OTHER  TEMPORARY  STACK  EXHAUSTED 

More  than  itAXTEM  temporary  storage  locations  are  needed  for  this 
statement 

unbalanced  parentheses 
unexpected  character  ’x* 

The  character  'X*  was  found  where  a  letter  was  expected 

UNRECOGNIZED  CONSTRUCT  ENCLOSED  BY  DECIMAL  °0lNTS 

UNRECOGNIZED  INSTRUCTION 

An  unrecognized  instruct  ion  to  the  CLUDGE  has  been  encountered 

UNRECOGNIZED  LIST  ITEM 

Improperly  constructed  list 

UNRECOGNIZED  STATEMENT 
Non-fatal 


UNRECOGNIZED  type 

An  unrecognized  type  occurs  in  a  *convert  instruction 

UNRECOGNIZED  type  BEGINNING  WITH  XXX 

A  type  na-ne  beginning  with  tbe  String  XXX  was  found  in  an  imp'ici 
type  declaration  statement 

WARNING  -  TYPE  OTHER  ARGUMENT  IN  I/O  LIST 
Non-fatal  (See  page  15) 
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Introduction 


This  report  I*  designed  to  serve  as  •  supplement  to: 

J.  M.  Yohe ,  "Best  Possible  Floating  Point  Arithmetic," 

The  University  of  Wisconsin,  Mathematics  Research  Center, 

U.S.  Army,  Technical  Summary  Report  March,  1970. 

A  description  is  given  here  of  the  implementation  of  "Best  Possible 
Floating  Point  Aritlwetic"  for  the  IBM  3&0/67  at  Washington  State  University, 
and  algoritlms  are  suppl ied  which  describe  the  Assembler  routines  used  in 
the  implementation.  A  complete  listing  of  these  routines  is  contained  in 
the  Interval  Arithmetic  Package. 


Indicator*  Used 


A.  Input  Ind icator: 


0PTI0N 


I*  the  neme  of  the  correction  option  Indicator.  Its 
values  and  results  are  as  folic***: 


0PTI0N 

l 

Sign  of  Result 

♦ 

1 

Least  Upper  Bound 

Augment  Truncate 

2 

Greatest  Lower  Bound 

Truncate  Augment 

3 

Truncate 

Truncate 

5 

Augment 

Augment 

U 

Round 

1 

Truncate  If  7th  hex  digit  of 
mantissa  of  result  is  ^8. 

Augment  otherwise. 

Interna  1 

Indicators: 

The  internal  indicators  used  are  essentially  the  sane  as  described 
in  MRC  Tech.  Report  #1051*. 


C.  Output  Indicator; 

FAULT  is  the  fault  indicator.  Its  values  and  their  meanings 
are: 


1 

E0 

exponent  overflew 

2 

INF 

inf  in i ty 

3 

EU 

exponent  underflow 

DZ 

division  by  zero 

Note:  0PTI0N  and  FAULT  are  the  first  and  second  words,  respectively, 
of  a  special  control  section  named  BPAIND  in  the  Assembler 


routines.  Hence  they  mey  be  accessed  from  Fortran  by  using  a 
labeled  ccmon  area  named  BPAIND.  or  from  Assembler  using  a 
dunvny  control  section.  OPTION  and  FAULT  are  both  treated  as 
Integer  quantities. 


originally  contain*  the  first  operand  in  addition,  the 
multiplicand  in  multiplication,  and  the  dividend  in  division. 


originally  contains  the  second  operand  in  addition,  the 
multiplier  in  multipi icetion.  and  the  divisor  in  division 


A  and  U  represent  single  precision  floating  point  numbers  as  they 
•re  represented  by  the  360,  l.e.,  the  high  order  bit  is  the  sign  bit, 
the  next  seven  bits  are  the  biased  exponent,  and  the  last  2^  bits  are 
the  normal ixed  mantissa. 


SIGNA 
SI  GNU  J 


these  contain  the  sign  bits  of  A  and  U  respectively 


r  x  pA 

y  these  contain  the  exponents  of  A  and  U  respectively 

EXPU 

J 

RESULT  contains  the  resulting  single  precision  f'oating  point 

number 


KANT  A  1^ 
KANTU  J 


these  contain  the  expanded  mantissas  of  A  and  U  respectively 


contains  the  multiplier  during  multiplication  and  the 
quotient  during  division 


HANTA,  KANTU,  and  V  are  double  words.-  In  the  algorithms,  MANTA,  and 
HANTA2  refer  to  the  leftmost  and  rightmost  words  of  KANTA,  and  KANTU, 
and  KANTU2  refer  to  the  leftmost  and  rightmost  words  of  KANTU. 


Below  are  given  the  algorithms  for  routines  BPAADO  (addition),  BPMUL 
(mult ipi teat  ion) ,  BPADIV  (division),  and  BPABND  (bounding). 

In  all  cases,  the  arrow  "« — **  indicates  replacement,  the  double 
arrow  indicates  mutual  replacement,  i.e.,  a  switching  of  values, 

and  a  subscript  of  16  indicates  a  hexidecimal  number.  Transfer  of  control 
is  indicated  by  the  next  step  to  be  executed,  without  the  words  "go  to." 
However,  if  transfer  of  control  is  to  a  different  routine,  the  words 


"go  to"  are  used. 


Addition  Algorithm  (BPAADD) 

1.  SC  ♦—  0; 

R1  ♦—  0; 

FAULT  ♦-  0; 

2.  If  |A  |  >  |u|,  A*-*  U; 

3.  If  A  4  0,  Step  5; 

U.  RESULT  ♦— U; 

Return; 

5.  If  A  t  0,  Step  7; 

6 .  A  ♦—  -A ; 


U  « - U 

SC  ♦-  1 


7.  itANTA  }♦—  0  ; 


*AKTA2*-  mantissa  of  A  with  00)6  as  rightmost  two  digits; 
EXPA «—  exponent  of  A; 

S  ICN'A  ♦-  sign  of  A; 


8.  MANTU ^  0 ; 

If  U  te  0,  Step  9; 

HANTUj  ♦-  one  's  complement  of  mantissa  of  U  with  FF^  as  rightmost 
two  digits; 

Step  10; 

9.  KANTU2«-  mantissa  of  U  with  00, 6  as  rightmost  two  digits; 

10.  EXPUe- exponent  of  U; 

S IGNU 4-  sign  of  U; 

11.  S COUNT  4-(EXPU  ‘  EXPA) x««; 

12.  If  SC0UNT  £  32.  Step  13; 


""W"' 


If  KANTA2  4  0,  Rl  ♦—  1; 

MANTAj  0; 

Step  14; 

13.  Shift  MANTA2  right  by  SC0UNT ; 

If  nonzero  digits  shifted  out,  Rl  «—l; 

14.  EX  PA  ♦—  EXPU; 

15.  MAMA  ♦-  MANTA  ♦  MANTU; 

16.  If  SICNU  4  0,  Step  20; 

17.  If  MANTA,  -  0,  go  to  BPABND; 

18.  Shift  MANTA  right  by  4; 

If  nonzero  digits  shifted  out,  Rle-1; 

>9.  EX  PA  ♦-  EX  PA  ♦  1; 

Co  to  BPABND; 

20.  if  MANTAj  4  0,  Step  22; 

21.  RESULT  *-0; 

Return; 

22.  If  Rl  *  1,  MANTA*-  MANTA  *  1; 

23.  MANTAj*—  one's  complement  of  MANTA^; 

SC  «—  1  -  SC; 

24.  If  first  hex  digit  of  MANTA^  4  0,  go  to  BPABND; 

25.  Shift  MANTA2  left  by  4; 

EXPA ♦—  EX  PA  -  1; 

Step  24; 

End 


n 


- r 

Add i t ion 


Multiplication  Al  oor 1 thm  (BPAMUL) 

1.  SC  4—  0; 

AI4-0; 

FAULT#-  0; 

2.  If  A  4  0,  Step  3; 

RESULT  4—  0; 

Return; 

J.  If  A  i  0,  Step  1«; 

A  4 - A; 

SC  4-1; 

U.  MANTAj  4-  0; 

MANTA  4-  mantissa  of  A  with  00j£  *s  rightmost  two  digits 
EXPA  ♦—  exponent  of  A; 

SICNA  4—  sign  of  A; 

5.  If  U  ¥  0,  Step  6; 

RESULT  4-  0; 

Return; 

6.  If  U  s  0,  Step  7; 

U  4 - U; 

SC  4-  1  -  SC; 

7.  KANTU (  4-  0; 

MANTU2  4-  mantissa  of  U  with  00^  as  rightmost  two  digits 
EXPU  4—  exponent  of  U; 

SICNU  4-  Sign  of  U; 

8.  V  4-  MANTA; 

MANTA  ♦-  0: 


9 


Hu 


t 

MC0UNT  ♦-  6; 

10.  Shift  V  r»9ht  by  *; 

11.  Shift  MANTA  right  by  A; 

If  nonzero  digit*  shifted  out,  M  «-  1; 

12.  RC0UNT  ♦-  15th  he*  digit  of  V; 

13.  If  AC0UNT  -  0,  Step  17; 

lit.  MANTA  «-  MANTA  ♦  MANTU; 

15.  RC0UNT  ♦-  RC0JNT  -  1; 

1b.  Step  13; 

17.  MC0UNT  ♦-  MC0UNT  -  1; 

If  MCJJNT  4  0.  Step  10; 

18.  If  MANTA,  -  0,  Step  20; 

19.  Shift  MANTA  right  by  A; 

If  nonzero  digit*  shifted  out,  Rl  ♦“  1; 
EXPA  «-  EXPA  ♦  1; 

20.  EXPA  «-  EXPA  ♦  EXPU  -  Al |fi; 

Co  to  BPABND; 

End 


9. 


If  MANTA  *  MANTU.  Step  11} 

MANTA  ♦-  MANTA  -  MANTU; 

v  ♦-  v  ♦  l; 

Step  9; 

||.  DC0UNT  OC0UNT  -  i; 

If  DC0JNT  -  0,  Step  13; 

12.  Shift  MANTA  left  by  4; 

Shift  V  left  by  U; 

Step  9; 

,3.  If  manta2  *  o,  AI  ♦-  1; 

1A.  MANTA  ♦-  V; 

15.  If  MANTA,  -  0,  Step  17; 

16.  Shift  MANTA  right  by  4; 

If  nonxero  digit*  shifted  out,  Rl  ♦-  1 
Add  1  to  EXPA ; 

17.  EXPA  ♦-  EX PA  - (EXPU  -  40, 6); 

Go  to  BPABND; 

End 


\ 


Bound ing  and  Round  1 ng  A 1 oo r I 1  hm  (BPABND) 

1.  If  EXPA  t  i*0|6.  Step  3; 

If  EXPA  £  7Fi6,  Step  6; 

2.  FAULT  ♦-  J; 

MANTA )  *-  0; 

HANTAj  *-  FFFFF FFO , 6; 

EXPA  ♦-  7F)6; 

Step  6; 

3.  If  EXPA  k  0,  Step  6; 

FAULT  «-  3; 

SCOUNT  ♦-  (-EXPA  ♦  5)*1*; 

If  SCOUNT  ^  32,  Step  U; 

If  MANTAj  ^  0 .  R I  a —  1 ; 

KANTAj  *-  0; 

Step  5; 

It,  Shift  MANTA  right  by  SCOUNT; 

If  nonzero  digits  shifted  out,  Rl  4—  1; 

5.  EXPA  ♦-  0; 

6.  If  rightmost  byte  of  MANTA  -  00)6  end  Rl  «  0,  Step 

7.  If  OPTION  •  3,  Step  Mi; 

If  OPTION  -  5,  Step  10; 

If  OPTION  4  It,  Step  8; 

If  rightmost  byte  of  MANTA  *  80,^,  Step  Mi; 

Step  10; 

8.  If  OPTION  -  J,  Step  9; 

If  SC  4  0,  Step  10; 


•3 


Step  >*♦; 

9.  If  SC  0,  Step  U*; 

10.  If  FAULT  4  3.  step  II; 

EX  PA  «-  0 
MANTA,  ♦-  0; 

MANTA2  ♦-  1 0000000, £  ) 

Step  ; 

!  I .  If  FAULT  4  1,  Step  12; 

FAULT  ♦-  2; 

Step  I1*; 

12.  MANTA  «—  MANTA  «•  1 00 | ^ ; 

If  MANTA,  -  0.  Step  1*4 ; 

13.  Shift  MANTA  right  by 

If  nonjero  digit*  ihifted  out,  Rl  *—  1; 

EXPA  EX  PA  ■»  1; 

If  EXPA  *  U0,6,  Step  U4; 

If  EXPA  *  7Fi&.  Step  2; 

IL.  If  SC  -  0,  Step  15; 

SICNA  4-  1  -  SICNA; 

15.  RESULT  «—  floating  point  number  whose  sign,  exponent, 

mentisse  ere  defined  by  SICNA,  EXPA,  MANTA 

Return; 


B round inc 


and 


End 
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c  r.  s  1  j  e 


e  shall  describe  a  method  for  computing  the 


r.ur.  value 


We  assure  the  number  of  points  in  X 


at  which  f  x 


are  a. wavs  correct 


sharr  these  bounds  car. 


lobal  c 


unction  ir. 


oerivatives  are  Known 


s  of  which  f  and  its  derivatives  are  defined  have  known 


ximations  w 


bounds  for  the  ar  curent  s  of  ir. teres 


?:nce  the  initial  box  can  be  chosen  as  larce  as  <e 


our  algorithm  actually  selves  the  unccnstrame 


rovided  it  is  known  that  the  solut; 


finite  region  (which  we  enclose  in  the  initial  box! 


There  is  a  corner  mi  sconcept ion  among  researcher*  in 
cr  t  it  i :  a  1 1  cr.  that  it  is  impossible  id  obtain  infallible  bounds 
on  x*  and  f*  computationally.  The  argument  is  that  we  can 
cr.lv  sample  :  \  and  a  feu  derivatives  of  fix  at  a  finite 
r.u-ber  of  points.  It  .s  possible  tc  interpolate  a  function 
having  the  necessarv  values  and  derivatives  values  at  these 
points  and  still  have  its  global  -inimum  at  ar.v  other  arbitrary 
point.  The  fallacv  of  this  argument  is  that  interval  ar.alvsis 
can  prcvide  bounds  or.  a  function  ever  an  entire  box:  that  is  cv 
continuum  of  points.  It  is  only  necessary  to  make  the  box 
sufficients  small  in  order  tc  make  the  bounds  arbitraril-  shar 
This  is  what  cur  algorithm,  does.  It  narrows  the  region  o: 
interest  until  the  bound  is  as  sharp  as  desired  subject  tc 
roundoff  restrictions  . 

In  a  previous  paper  ;3],  »e  gave  a  method  of  this  tvpe  for 
the  one  *  d  i  ~en  s  i  or.a  I  case.  The  method  never  failed  tc  converge 
provided  f  x  and  f  x  had  on!"  a  finite  number  of  isolated 
oercs.  Our  method  for  the  n-di-er.s  tonal  proble-  appears  tc 
alwavs  converge  also;  but  i.e  have  net  vet  attempted  to  prcv< 
•sher.  it  does  converge,  there  is  never  a  question  that  x*  and 
satisfy  the  computed  bounds. 

Recently,  R .  E.  Moore  [ID]  published  a 
the  range  of  a  rational  function  of  r.  variables  over  a  bounded 
region.  (See  also  [14].)  Although  he  does  not  note  the  fact, 
his  method  will  serve  to  bound  the  global  minimum  value  f *  of 
a  rational  function.  However,  cur  algorithm  is  more  efficient. 
Moreover,  it  is  designed  to  bound  x*  as  well  as  f* . 


emr  t  e d  t  : 

prove  it 

icr.  that 

x*  and  f* 

thod  for 

computing 

es  over  a 

bounded 

he  suggest  the  reader  read  the  previous  paper  [5]  before  the 
current  one.  The  one-dim.ens  ional  :3>e  therein  serves  as  ar.  eas.er 
introduction.  However,  the  current  paper  is  esser.tiallv  self 
contained.  It  would  be  better  if  the  reader  had  some  familiarity 
with  the  rudiments  of  interval  analysis  such  as  car.  be  found  :n 
the  first  three  chapters  of  [9).  However,  we  shall  review  some 
of  :ts  relevant  properties. 

Our  method  will  find  the  global  -.ini mum  or  minima''.  Because 
of  computer  limitations  of  accuracy,  it  r.a"  also  find  near-global 
minima  such  that  rounding  errors  prevent  determination  of  which 
is  the  true  minimum.  However,  if  the  ter-ir.atirr  criteria  are 
sufficiently  stringent,  our  algorithm  .ill  alwavs  eliminate  a 
local  minimum  whose  value  is  substantial ly  larger  than  f*  . 

,ur  algorithm  is  composed  of  four  separate  parts.  'ne  part 
uses  an  interval  version  of  Newton’s  method  to  find  «t*::?nar- 

f  \ 

points,  A  second  pa.t  eliminates  points  of  N  where  f  is 
greater  than  the  smallest  currently  known  value  T. 

A  third  part  of  our  algorithm  tests  whether  :  is  monotonic 
in  a  sub-box  X  of  X  '  T  (  *•-  '’■*’**“  r'nT*  r'r  a  1 1  n*  v 


If  so,  we  delete  part  or  ail  of  X 


depending  on  whether  X  contains  boundary  points  of  X  . 

A  fourth  part  checks  for  convexity  of  f  in  a  sub-box  X  of  X 


If  f  is  not  convex  anvwhere  in  X,  there  cannot  be  a  stationarv 
..point  of  f  in  X. 


The  first  part  cf  the  algorithm,  if  used  alone,  would  find 
all  stationary  points  in  X  .  The  second  part  serves  to 
eliminate  stationary  points  where  f  >  f*  .  Usually  thev  are 
eliminated  before  thev  are  found  with  any  great  accuracy.  Hence 
computational  effort  is  net  wasted  using  the  first  part  to 


ne  se. 


accurate!  find  ar.  unwanted  stationary  point, 
also  serves  to  eliminate  boundarv  points  of  X  and  to  find  a 
global  minimum  if  it  occurs  on  the  boundary.  The  second  part  cf 
the  algcrith-  used  3l:r.e  would  find  the  global  -inimum  or  minima 
but  its  asymptotic  convergence  is  relatively  slow  compared  to  that 
of  the  Newton  method.  Hence  the  latter  is  used  also.  The  third 
and  fourth  parts  of  the  algorithm  merely  improve  convergence. 

2  .  I  r.terva  1  analysis  . 

The  tool  which  allows  us  to  be  certain  we  have  bounded  the 
global  mm imu-  is  interval  analysis.  We  bound  rounding  errors  bv 
using  interval  arithmetic.  "ere  importantly,  however,  we  use 
interval  analysis  to  bcum.d  the  range  c:  3  function  over  a  hex. 

Let  g(x  be  a  rational  function  cf  n  variables 
On  a  computer,  we  can  evaluate  g(x  for  a  given  x  by 
performing  a  sequon^e  of  arithmetic  operations  involving  only 
addition,  subtraction,  multiplication,  and  division. 

Let  X.  (i  •  l,...,n)  be  closed  intervals.  If  we  use  X. 
in  place  of  x  and  perform  the  same  sequence  of  operations  using 
'interval  arithmetic  (s«e[9j)  rather  than  erdinarv  real  arithmetic, 
we  obtain  a  closed  interval  g  X)  containing  the  range 


(g(x)  :  x.  *  Xi  (i  •  1 . n)> 


o 


o:  g.x)  over  the  box  X  .  This  result  will  not  re  sharp,  m 

general,  but  if  outward  rounding  (see;.-]  is  used,  then  g(X 
will  always  contain  the  range.  The  lack  cf  sharpness  results 
from  other  causes  besides  roundoff.  With  exact  interval 
arithmetic,  the  lack  cf  sharpness  disappears  as  the  widths  of 
the  intervals  decrease  to  cero. 

If  g(x)  is  not  rational,  we  assume  an  algorithm  is  known 
for  computing  an  interval  g  X-'  containing  the  range  of  g(x 
for  x  -  X  .  Methods  for  deriving  such  algorithms  are  discussed 
in  (9]). 

3 .  Taylor's  Theorem. 

We  shall  use  interval  analysis  in  con; unction  with  Taler’s 
theorem  ir.  two  wa--s.  First,  we  exrand  f  as 


(3.1  • 


f  y)  •  f  x;  *  (y*l  1  5  x  •  i(y-x)T  H  x.v.?}  (v-x^ 


where  g'x  is  the  gradient  of  f  x)  and  has  components 
g,  \  *  \  :-x.  .  The  quant itv  H  x,v,-;  is  tr.e  Hessian,  matrix 

to  be  defined  presently.  For  reasons  related  to  the  use  of  interval 
ar.alvsis,  we  shall  express  it  as  a  lower  triangular  matrix  instead 
of  a  symmetric  matrix  so  that  there  are  fewer  terms  in  the  quadratic 


form  involving  H(x,y,£)  . 

Ke  define  the  element  in  position  (i,j)  of  Hfx.v,;)  as 


/e-f/ax."  for  i  •  i  i-l,***,n)  , 

•  |  2J"f/5xi?x  .  for  j<  i  (i  ■  1,  •  ••  ,n;  j  ■  1,  ••  •  ,i-l) 


for  j  •  i  i  -  1,  ••  •  ,n>  , 


:3.2)  h.. 


0  otherwise. 


The  arguments  of  hi-  depend  on  i  and  i  .  If  we  expand  f 
sequentially  in  one  of  its  variables  at  a  time,  we  can  obtain  the 


Then 


Assume  x.  f  X.  an 


the  arguments  of  H 


for  each  i  •  1 


For  general  n 


Other  arrange-ents  of 


point  in 


or  ar.v  pc  in 


In  the  secuel,  we  sha 


ter  notation  anc  jse 


a?  nanv  argu’-e 


nterva,  quantities 


of  h':>  as  possible.  The  standard  Tavlor  expansion  would  have 


ntervals  for  all  arguments  of  all  elements  o 


A  more  general  approach  o 


of  expansion  was  introduced  in 


he  other  Tavlor  expansion  we  shall  want  is  of  the  gradien 


n)  of  g  can  be  expanded  as 


Each  element 


where 


atrix 


How eve 


e  are  exran 


H  is  lower 


the  elements  o 


elenen 


he  same  arguments  a? 
e  need  only  calculate 


have 


follows  easily 


ate  value  cf  the 


r.um 


eletes  sub -boxes  of  X 


A  test  for  convexity. 


S . 

As  our  algorithm  proceeds,  we  dynamically  subdivide  X  ' 
into  sub-boxes.  Let  X  denote  such  a  sub-box.  •'>  evaluate 

*  (Xl  •  *'*  »xn*J  for  i  ■  1.  •••  ,n  ,  where  h^  is  the  diagonal 
element  of  the  Hessian.  Note  that  every  argument  of  h , ,  is  an 
interval  and  hence  the  resulting  interval  contains  the  value  cf 
h ,  i  ( x  ’  for  every  x*X  .  That  is,  if  lUj.vJ  denotes  the  compu¬ 
ted  interval  h^fX,,  •••  ,Xn)  ,  then 


fcr  all  x  *  X  . 


Suppose  we  find  v.  <  0  for  some  value  of  i  .  Then 
h ;  j  ( x )  <  0  for  every  x<X  .  Hence  there  is  no  point  in  X 
which  the  real  non- interval)  Hessian  is  positive  definite. 


at 

Hence 


f  is  not  convex  and  cannot  have  a  minimum,  which  is  a  stationary 
point  in  X  .  Hence  we  can  delete  all  of  X  except  for  any  boun¬ 
dary  points  of  X^0)  which  might  lie  in  X  . 

When  we  evaluate  h^(Xj,  •••  ,Xn)  ,  we  may  find  that  the  left 
endpoint  uiO  for  all  i  •  1.  •••  ,n  .  when  this  occurs,  we 
know  from  inclusion  monotonicity  (see  [9])  that  we  will  find  each 


>  o  for  anv  sub-box  of  X  .  Hence  we  could  save  some  computa¬ 
tional  effort  by  noting  when  a  box  is  a  sub-box  of  one  for  which 
u  >  C  for  all  i  •  1,  •••  ,n.  *e  would  skip  this  test  for  such  a  box. 

Note  that  an  element  h  ^  with  arguments  (Xj,  **•  ,Xn)  is 
not  obtained  when  we  compute  H  X)  since  the  diagonal  elements  of 
H (X)  have  arguments  different  from  (X..  •••  ,Xn)  except  for  the 
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ele-ent  in  position  (n,n).  Hence  our  test  for  convexity 
requires  recalculation  of  the  diagonal  of  the  Hessian. 

6  •  The  interval  Newton  method  . 

For  each  sub-box  X  of  X  '  that  our  algorithm  generates, 
we  car.  appiv  ar.  interval  \ewton  method  to  the  gradient  g.  Such 
methods  seek  the  ceros  of  g  and  hence  the  stationary  points  of 
f.  Such  a  method  procudes  from  X  a  new  box  or  boxes  N^X  .  Any 
points  in  X  not  in  N,X  cannot  contain  a  cere  of  g  and  car.  be 
discarded  unless  they  are  boundary  points  of  X  '  . 

These  methods,  in  effect,  solve  ;'3.5)  for  points  v  where 
g  >•)  •  0.  The  first  such  method  was  derived  by  Moore  (9). 

Variants  of  Moore's  -ethod  can  be  found  in  [3,6,8,13).  The  most 
efficient  variant  can  be  found  in  [6].  Krawccyk's  method  [S’ 
is  a  suitable  alternative  to  the  method  in  (6).  Discussions 
of  Krawccyk’s  method  can  be  found  in  [11]  and  (11). 

U'e  now  give  a  brief  svnopsis  of  the  method  in  (6).  he  wish 
to  solve  the  set  of  equations 

g(x)  *  J(CHy*>  •  o 

for  the  set  of  points  y  obtained  by  letting  C  range  over  X  . 
he  shall  find  a  subset  V  of  X  containing  this  set. 

Let  J.  be  the  matrix  whose  element  in  position  (i,i'  is 
the  midpoint  of  the  corresponding  interval  element  :tjm  of  th. 
Jacobian  J(X)  .  Let  B  be  an  approximate  inverse  of  J£  .  As 
pointed  out  in  [3],  a  useful  first  step  in  solving  for  Y  is  to 
multiply  (6.1)  by  B  giving 
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(6.2) 


Bgfx)  *  BJCi)(y*x)  -  0  . 

Note  that  the  product  BJ(£)  approx i: nates  the  identity  matrix. 
However  it  may  be  a  very  poor  approximation  when  X  is  a  large 


■6.3) 


We ’*olve"  (6.2  by  a  process  similar  to  a  single  sweep  cf  the 
Gauss  -  Seidel  method.  Write 

BJ(X)  •  l  *  D  *  U 

where  L,  D.  and  U  are  the  lower  triangular,  diagonal,  and 
upper  triangular  part  of  BJ(X  ,  respectively.  .he  inter\al  a...x 

0  *  ■  diag[l  G>.,  •••  ,1/-^ 

contains  the  inverse  of  every  matrix  in  D  .  -he  box  i  sol\.ng 
(6.2)  is  obtained  as 

Y  -  x  *  D*1 (Bg(x)  *  LC^*X)  *  U(X*x) )  . 

When  obtaining  the  component  of  )  ,  the  components  N.j,  ’^j-1 

appearing  in  the  right  member  of  this  equation  have  already  been 

obtained . 

This  formulation  presupposes  that  the  intervals  fi  •  1,  •  •  •  .n 

•■do  not  contain  zero.  When  X  is  a  small  box,  BJ  ( X  is  close.- 
approximated  by  the  identity  matrix  and  hence  D  is  also.  However, 
for  X  large,  it  is  possible  to  have  0  f  Dii  for  one  or  more 
values  of  i  .  This  case  is  easily  handled.  He  simply  use  ar. 
extended  interval  arithmetic  which  allows  division  by  an  interval 
containing  zero.  (See  (6]  for  details.) 


Note  that  we  cannot  allow  the  Newton  procedure  to  delete 
boundary  points  of  X  '  since  the  global  minimum  need  not  be  a 
stationary  point  if  it  occurs  on  the  boundary.  K'e  discuss  this 
point  further  in  Section  10. 

If  we  were  to  use  this  Newton  method  only,  w*  would  :r. 
general  find  stationary  points  of  f  which  were  not  minima. 
Vcreover ,  we  would  find  local  minima  which  were  not  global 
minima.  To  avoid  this,  we  use  an  additional  procedure  to 
delete  points  where  f  exceeds  the  smallest  known  value  7.  This 
procedure  is  described  in  the  next  section. 

In  some  applications,  it  nav  be  desirable  to  find  all  the 
stationary  points  of  f  ir.  a  given  box.  This  can  be  dene  using 
the  Newton  method  alone  or  in  conjunction  with  the  mcnotonicit- 
check  of  Section  9.  If,  in  addition,  the  convexity  check  cf 
Section  5  were  used,  all  stationary  points  except  maximum  would 
be  found. 


MMmmu 
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(M) 


7 .  Bounding  f  . 

*>  now  consider  how  to  delete  points  y  <X  where  we  know 
f  v)  >  f  and  hence  where  fly)  is  not  a  global  minimum.  -e 
retain  the  complementary  set  which  is  a  sub-box  or  sub-bcxes 

Y c  X  wherein  f(y)  "ay  be  5  f  . 

As  pointed  out  in  [ S] .  if  we  only  wish  to  bound  f*  and  net 

x*  ,  we  car.  delete  points  where 

f(y)  >  ?  - 

for  some  e,  >  0  .  he  can  allow  Cj  to  be  nonzero  only  if  we  do 
not  need  to  know  the  pcmt(s)  x*  at  which  f  is  globally  mini¬ 
mum  . 

we  want  to  retain  points  where  .‘.1  is  not  satisfied.  Fro.. 
(3.1),  this  is  the  case  for  points  y  if 

f(x)  *  (y-x)T  g(x)  *  i(y-x)  H  ;  y-x  -  f  *  ci 

because  the  left  member  equals  fly)  •  Tenote 

E  •  f  *  f(x)  *  tj  . 

Then 


(T.2) 


yTg(x)  ♦  yyT  H(Oy  s  E 

where  y  •  y-x  .  We  shall  use  this  relation  to  reduce  X  m  one 
dimension  at  a  time  to  yield  the  sub-box(ts)  Y  resulting  from 
deleting  points  where  f(y)  >  ?  -  c1  • 


We  shall  illustrate  the  process  for  the  case  n  ■  *  .  The 
higher  dimensional  case  follows  in  the  same  wav.  For  n  ■  2  , 

(".21  becomes 

*  Y;g:(x)  ♦i|yf  hu(0  *  yjVVUj (O  ♦  Vh22(£-I  5  E  • 

'.Ve  first  wish  to  reduce  X  in  the  x ,  -  di  rec  t  ion  .  Thus  we 
solve  this  relation  for  acceptable  values  of  y ^  .  After  collecting 
terms  in  v,  ,  we  replace  y,  by  X,  .  In  the  higher  dimensional 
case  we  would  also  replace  y  by  for  all  i  •  3,  •••  ,n  .  he 

also  replace  i  by  X  (since  C  <  X)  .  We  obtain 

y.,lg1  x  4X,h,.(X)]  *  iy^ hxl (X'  *  XjgjCx)  *  ii;*h22W  *  E  5  0 
wnere  X,  •  X,  *  x„  . 

We  solve  this  quadrat::  for  the  interval  or  intervals  of  points 
V as  described  below.  Call  the  resulting  set  .  Since  we  are 

only  interested  in  points  with  y ^  <  Xj  ,  we  compute  the  desired 


set  Y  j  as  Y  ^ 


y  r  ? 

*1  . 


For  the  sake  of  argument,  suppose  Y ^  is  a  single  interval. 

We  can  then  try  to  reduce  X,  the  same  wav  we  (hopefully)  reduced 

m 

X,  to  get  Yj  .  We  again  rewrite  (*.3).  This  time  we  replace 

by  Yj  and  i.as  before)  C  by  X  .  We  could  obtain  better  results 
by  replacing  ^  bv  Yj  rather  than  Xj  but  this  would  require 
re-evaluation  of  the  elements  of  H  .  We  obtain 

?;{g;(x)4?ih;i(X)]  *  7y22h22(X)  *  VjgjCx)  *  ^f'hjjCX)  -  E  *  0 


15 


If  the  solution  set  V,  is  strictly  contained  in  X;  ,  we 

could  replace  X,  by  Y,  in  ('.4)  and  solve  for  a  new  Yx  .  *e 
•»  • 

have  not  tried  tc  do  this  in  practice.  Instead,  we  start  over 
with  the  box  Y  in  place  of  X  as  soon  as  we  have  tried  to  reduce 


each  Xx  to  Yi  (i  ■  1 . 


,n?  .  Note  th;s  r.ear.s  we  re-evaluate 


h;x} 


v»e  now  consider  how  to  solve  the  quadratic  equation  •  .  o. 

('.5).  These  have  the  general  form 


C.6) 


A  *  St  *  Ct*  5  0 


where  A  .  B  ,  and  C  are  intervals  and  we  seek  values  of  t 
fying  this  inequality. 


satis 


Denote  C  •  (c^c,)  and  let  c  be  an  arbitrary  point  in  C  . 


Similarly,  let  a«A  and  b«B  be  arbitrary.  Suppose  t  is  such 


that  (,“.e)  is  violated;  that  is  Q(t)  >  0  ,  where 


Q(t)  •  a  ♦  bt  *  ct”  . 


If  this  is  true  for  c  ■  Cj  ,  then  it  is  true  for  all  c  <  C  . 


Hence  if  we  wish  to  find  the  complementary  values  of  t  where  (  .*) 
might  hold  we  need  only  consider 


(7.7) 


A  ♦  Bt  *  Cjt‘  «  0  . 


If  Cj  •  0  ,  this  relation  is  linear  and  the  solution  set 

is  as  follows:  Denote  A  •  and  B  "  ^bl  ,b2^  *  Thcr‘  The 

set  of  solution  points  t  is 


f-aj/b;,-]  if  «  0  ,  b,  <  0  , 

r-aj/bj.-J  if  aj  >  0  ,  bj  <  0  ,  b,  5  0  , 

['*•*)  if  »j  i  0  ,  b,  :  0  <  b,  , 

I-.-aj/bjJ  .  [-a1/b1#-]  if  a1  >  0  ,  bj  <  0  <  b,  , 

t— .-.j/bj)  if  ,J  <  o  ,  bj  »  0  , 

i^.-a^/b,]  if  ai  >  0  ,  b^  »  0  ,  b-  :>  0  , 

enptv  set  if  a^  >  0  ,  b.  “  b-.  ■  0  . 

Recall  that  we  will  intersect  T  with  X,  for  some  value  of  i  . 

Thus  although  T  may  be  unbounded,  the  intersection  is  bounded. 

If  c,  *  0,  the  quadratic  C".6)  may  have  no  solution  or  it 
may  have  a  solution  set  T  composed  of  either  one  or  two  ir.terva 
Tn  the  latter  case,  the  intervals  mav  be  semi  -  inf inite .  However 
after  intersecting  T  with  X,,  the  result  is  finite. 

Denote 

Qx( t)  -  a  *  bt  *  Cjt” 

where  a  (  A,  b  <  B.  and  cx  is  the  left  endpoint  of  C.  We  sha 

delete  ooints  t  where  Q^t)  >  0  for  all  a  f  A  and  b  <  B. 
Thus  we  retain  a  set  T  of  points  where  Qx(t)  ~  P«  as  desired. 

But  we  also  retain  (in  T)  points  where,  for  fixed  t,  QjU  >  0 

for  some  a  i  A  and  b  (  B  and  Qx(t)  i  0  for  other  a  <  A  and 

b  t  B.  This  same  criterion  was  used  to  obtain  T  when  ■  '  . 
This  assures  that  we  shall  always  retain  points  in  x i  where 
f(x)  is  a  minimum. 


Denote 


q^t)  - 


q,U) 


b;t 

*  Clt2 

if 

t 

< 

0, 

V 

4  clt 

if 

t 

> 

0 

V 

*  Cit2 

if 

t 

< 

0, 

7 

b2t 

*  clr‘ 

if 

t 

> 

0. 

Then  we  can  write  the  interval  quadratic  as 

Qx  C t )  ■  [ara:l  -  (bj.bjlt  *  cxt* 

•  [q1ft),  q 2 ( t ) ]  . 

Thus  for  any  finite  t,  q^^  ( r )  is  a  lower  bound  for  C,  (t)  and 
q,(t)  is  an  upper  bound  for  Q,ft)  for  any  a  <  A  and  any  b  <  B. 

For  a  given  value  of  t,  if  OjCO  >  0,  then  Qj  ft)  >  0 
for  all  a  *  A  and  b  i  B.  Hence  we  need  only  to  solve  the  real 
quadratic  equation  q.(t  •  0  in  order  to  determine  intervals 
wherein,  without  question,  Q, (t)  >  0.  This  is  a  straightforward 
problem . 

The  function  q.(f  is  continuous  but  its  derivative  is 
discontinuous  at  the  origin  when  *  b,  which  will  generally 
be  the  case  in  practice.  Hence  we  must  consider  the  cases 
t  ;  0  and  t  *  0  separately. 

If  c 1  >  0,  the  curve  q1(t)  is  convex  for  t  <  0  and 
convex  for  t  >  0.  Consider  the  case  t  «  0.  If  q ^ ( t )  has  real 
roots,  then  Qj(t)  >  0  outside  these  roots,  provided  t  :  0. 

Hence,  we  retain  the  interval  between  these  roots.  Ke  need  only 
examine  the  discriminant  of  c.(t)  to  determine  whether  the  roots 
are  real  or  not.  Hence  it  is  a  simple  procedure  to  determine  which 
oart  (if  anv)  of  the  half  line  t  i  0  can  be  deleted.  The  same 
procedure  can  be  used  for  t  ;  0. 


For  c1  <  0,  qj(t)  is  concave  for  t  s  0  and  for  t  >  0. 

In  this  case  we  can  delete  the  interval  (if  any)  between  the  roots 
of  c ^ ( t )  in  each  half  line.  The  set  T  is  the  complement  of 
this  interval.  It  is  composed  of  two  semi  -  inf ini te  intervals. 

In  determining  T  for  either  the  case  c1  <  0  or  in  the 
case  c,  >  0,  it  is  necessary  to  know  whether  the  discriminant 


is  non-negative  or  not.  Denote 


when  t 


These  are  the  discriminants  o 


respect ively . 

When  we  compute  A.  or  A , ,  we  shall  make  rounding  errors 
Thus  we  should  compute  them  using  interval  arithmetic  to  boun 
these  errors.  When  comnuting  A.  •  (i  ■  !»»)•  suppose  we  obt 


the  interval 


We  use  the  appropriate  endpoint  of  Aj  or  A-  to  determine 
which  assures  that  we  never  delete  a  point  t  where  Qj't  couaC 
be  non - po s i 1 1 ve  .  Thus  we  use  the  endpoint  of  A^  or  vhi.h 

yields  the  larger  set  T  . 

When  we  commute  the  roots  of  q^(ti,  w«*  s h a  1 1  make  rounding 
errors.  Hence  we  compute  them  using  interval  arithmetic  and  again 
use  the  endpoints  which  yield  the  larger  set  T  to  assure  we  do 
not  delete  a  point  in  X.  where  f  is  a  minimum. 


denote 


As  is  well  known,  the  rounding 


Note  that  R 


in  the  form  R.  rather  than  in 


error  is  less  if  we  compute  a  too 
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the  form  S’  when  <  0.  The  converse  is  true  when  bi  >  0. 
Similarly,  the  rounding  error  is  less  when  using  R^  rather  than 
S*  when  >  0.  Hence  we  compute  the  roots  of  q  ^  C  t )  as  Ri 
and  S*  when  b ,  <  0  and  as  Rj  and  S:  when  b^  >  0. 

Note  that  computing  R^  or  involves  taking  the  square 

root  of  the  interval  C*  .  In  exact  arithmetic  this  would  be  the 
real  quantity  We  would  never  be  computinc  roots  of  qj(t) 

when  _ i  was  negative.  Hence  if  we  find  that  the  computed  result 
C 1  contains  rerc,  we  can  replace  it  by  its  non-negative  part. 

Thus  we  will  never  try  to  take  the  square  root  of  an  interval 
containing  negative  numbers. 

Giver,  any  interval  I,  let  IL  and  I*  denote  its  left  and 
right  endpoint,  respectively.  We  use  this  notation  below.  Ising 
the  above  prescriptions  on  how  to  compute  the  set  T  ,  we  obtain 
the  following  results: 


For 

bl  * ( 

[>  and 

Cj  >  0. 

j 

:  (the 

empty  set 

if 

R 

“2  ‘ 

:  0. 

(?. 

9)  T  ■ 

(s'2)r3 

if 

a1  » 

0  and 

“2  5 

1 

{ 

(s;)ri 

if 

al  S 

0. 

For 

bj  :  C 

and  c1 

>  0. 

j 

f  * 

if 

0. 

v- 

9)  T  - 

K>R] 

if 

al  > 

0  and 

t> 

— *  30 

IV 

o 

• 

j 

\ 

(Rj)R3 

if 

a1  5 

0. 

Note  that  if  c1  >  0,  then  l2  car  be  negative  only  if  a]  >  0.  uence 
the  condition  l2  <  0  implies  a^  >  0.  This,  and  similar  cases,  has  been  used 
to  shorten  the  conditional  statements  in  the  above  expressions  for  T  . 

We  have  seen  that  the  solution  of  the  quadratic  inequalities  such  as 
(7.4)  or  (7.5)  can  be  an  interval  or  two  semi- infinite  intervals,  say 

and  z/2)  .  The  desired  solution  set  ^  is  obtained  by  intersection 
with  X1  .  In  the  former  case,  Y^  «  X^  1  Z^  can  be  empty  or  a  single 
interval.  In  the  latter  case. 


can  be  empty,  a  single  interval,  or  two  intervals.  Ke  now  consi¬ 
der  the  logistics  of  handling  these  cases. 

The  quadratic  inequality  to  be  solved  for  Z,  will  have 

quadratic  term  -i-y^h^CX)  so  the  interval  C  in  (“.6)  is 

i  h%i(X'  and  the  left  endpoint  is  c.^  -  [i  h^(X)]*”  .  If 

c^'1'  >  0  the  solution  set  is  a  single  interval.  But  if  c^  <  0, 

it  is  two  semi- infinite  intervals  and  it  may  be  that  will  be 

two  intervals.  This  would  complicate  the  process  of  finding 

Y.  ,  ,  •  •  •  ,Y  .  Thus  we  proceed  as  follows. 
i*i  n 

Let  I.  denote  the  set  of  indices  i  for  which  c^1  ?  0 
and  I,  denote  the  set  of  indices  i  for  which  <  0  .  He 

first  find  Y^  for  each  i  *  Ij  .  Ke  then  begin  to  find  Y^  for 

i  «  I,  .  Let  j  €  I,  be  such  that  Y.  is  composed  of  two  inter- 
vals,  say  Y.  1  and  Y.  '  .  Then  X,  is  the  smallest  interval 

containing  both  Yi  *  and  Y ,  .  When  finding  Y^  for  the 

remaining  values  of  i  ,  we  use  X 1  in  place  of  Y ^  . 

After  finding  all  Y^  for  i  ■  1,  •••  ,n  ,  we  wish  to  use  the 
fact  that  we  can  delete  the  interval,  say  Y,c  ,  between  Y. 

and  Y.^  .  We  would  like  to  do  this  for  all  the  values  of  j 

for  which  Y ^  was  two  intervals.  However,  it  could  be  that  this 
occurred  for  1  the  indices  j  •  1,  ,n  .  After  deleting  the 
interior  interval  Y.c  in  each  dimension,  the  resulting  set  would 
be  composed  of  2n  boxes.  For  large  n  ,  this  is  too  many  boxes 
to  handle  separately.  Hence  we  delete  only  a  few  (one,  two,  or 
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three  of  the  largest  of  the  intervals  Y.c  .  We  then  process 
each  of  the  new  boxes  separately. 

We  would  like  to  prevent  the  generation  of  long,  narrow 
boxes.  Thus  a  good  choice  of  which  Y^  to  delete  is  the  cr.e  s 
correspond ir.g  to  the  conpcnent  for  which  the  smallest  interval 
containing  both  Y.  1  and  Y i v  *  is  largest.  However,  we  have 
chosen  to  delete  the  largest  interval  Y,w  . 

Let  us  call  the  process  we  have  described  in  this  section 
the  We  car.  combine  the  quadratic  method  with 

the  Newtcr.  method.  It  is  desirable  tc  do  this  35  we  now  explain. 

If  the  left  er.dpomt  of  H,  .  X  is  negative,  then  the  quad- 

X  * 

i  1  ^  I 

ratio  method  car.  give  rise  to  two  new  intervals  Y.  ‘  ar.d  Y. 
in  place  of  X.  .  When  trying  to  improve  X.^,  (say),  it  is 
impractical  tc  use  Y  ■  '  *  and  Y,1*  separately  and  we  use  X.  , 

instead.  Thus  the  improvement  of  X.  is  of  nc  help  when  trying 
to  improve  X .  ^  »  etc.  Similarly,  when  applying  the  Newton 
step,  if  Jii(X)  contains  tero  as  an  interior  point,  we  can 
obtain  two  subintervals  in  place  of  X i  .  Again,  we  canr.ct  con¬ 
veniently  use  this  fact  in  the  remaining  part  of  the  Newton  step. 

We  wrould  like  tc  dc  these  steps  first  whi~\  are  c:  help  in 
subsequent  steps.  Hence  the  following  sequence  is  suggested. 

First  try  to  improve  X ^  by  the  quadratic  method  for  each  value 
of  i  •  1,  •••  ,n  for  which  the  left  endpoint  of  H^(X)  i? 
positive.  Then  apply  the  interval  Newton  method  to  the  (old  cr 


( 


new)  components  for  which  (M3  J^(X)(i  ■  1,  •••  ,n)  .  Next  use 
the  quadratic  method  for  these  compcr.er.ts  for  which  the  left  end¬ 
point  of  H^(X  1S  non-positive.  Finally,  complete  the  Newton 
step  for  those  components  with  CM  B  J^.(X)  . 

.At  each  stage  of  either  method,  when  trying  to  improve  the 
i-th  component  of  the  box,  we  use  the  currer.tlv  best  interval  fer 
the  other  components.  This  may  be  the  smallest  interval  contain¬ 
ing  two  disjoint  intervals  in  some  cases.  In  fact  it  would  be 
possible  for  the  quadratic  method  and  the  Newton  method  to  each 
delete  disjoint  sub  -  interval  s  for  a  giver,  component .  This  would 
give  rise  to  three  sub- intervals  to  be  retained.  However,  it 
seems  better  to  simplify  this  case  and  only  delete  the  larger  of 
the  two  sub- intervals . 

When  both  methods  are  completed,  we  may  have  several  compo¬ 
nents  divided  into  two  sub- intervals .  If  so,  we  find  the  one  for 
which  the  largest  interior  sub-interval  has  been  deleted.  We 
replace  all  the  others  by  the  smallest  sub-interval  containing 
the  two  disjoint  parts.  We  then  divide  the  remaining  part  of 
the  current  box  into  two  sub-boxes  bv  deleting  the  sub-interval 
for  the  component  in  question.  Ke  could  do  this  for  more  than 
one  component,  but  each  deletion  would  double  the  number  of  boxes. 
It  seems  better  to  keep  the  number  of  boxes  small. 
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5.  Choice  of  e,  • 

Suppose  we  want  to  bound  the  value  f*  of  the  global  mimaun 
to  within  a  tolerance  tj  but  we  dc  not  care  where  :  takes  cr. 
this  value.  Then,  as  pointed  out  in  Section  ~,  we  can  delete 
points  y  where 

f(y)  >  f  -  tj  • 


Once  we  have  found  a  point  x  where  £(x) 
f  -  f*  «  e,  ,  our  algorithm  will  eventual 
if  we  use  (8.1).  However,  x  r.av  be  far 

f  A  \ 

f  is  globally  minimal.  When  all  of  X 


*  f  is  such  that 

(O', 

ly  delete  all  cf  X 
fror  the  point  x*  where 
i«  deleted,  we  will 


know  that 


f  -  tj  i  f*  !  f  . 

Choosing  c,  >  0  will  speed  up  our  algorithm.  However,  if 
we  wish  tc  obtain  good  bounds  on  x*  ,  we  must  choose  Ej  "  •  * 

We  then  terminate  our  algorithm  when  the  remaining  set  or  points 
is  sufficiently  small.  See  Section  13  for  a  termination  procedure. 


9.  Monotonicitv . 


Another  step  in  our  algorithm  makes  use  of  the  monotonicitv 
of  f  .  Suppose,  for  example,  the  i-th  component  g.( x'  of  the 
gradient  is  non-negative  for  all  x  1  X  .  Then  the  smallest  value 
of  fx)  for  x*X  must  occur  for  x.  equal  to  the  left  end¬ 
point  of  X^^  . 

To  make  use  of  a  fact  such  as  this,  we  evaluate  g,(X^,  **•  ,Xn) 
The  resulting  interval,  which  we  denote  by  ,  contains 

gi (x)  for  all  x<X  .  Denote  X.  ■  [xi“,x,^]  .  If  >  0  , 
f^x)  is  smallest  (in  X  )  for  x}  •  x,*-  .  Hence  we  can  delete 
alio:  X  except  the  points  with  x.  m  X1L  •  If  >  0  ,  fix' 
cannot  have  a  stationary  point  in  X  .  Hence  we  can  delete  ill 
of  X  unless  the  boundary  at  x,  *  x . ^  contains  boundary  points 
of  the  initial  box  Xv  ^  .  Similar  results  occur  if  «.  ■  t  0  or 
if  wi  <  0  . 

We  evaluate  g^Xj,  *’*  ,Xn)  for  i  -  1,  •••  ,n  and  reduce 

the  dimensionality  of  X  for  any  value  of  i  for  which  o^  i  0 

or  ^  «  0  .  Of  course,  we  delete  all  of  X  ,  if  possible. 

It  is  possible  that  we  can  reduce  X  in  every  dimension  in 

this  wav.  If  so,  only  a  single  point,  say  x  ,  remains.  In  this 

case,  we  evaluate  f(x)  .  If  f(x)  >  f  ,  we  can  eliminate  x 
and  hence  all  of  X  is  deleted  by  the  process.  If  f(x)  s  f  , 
we  reset  7  equal  to  f(x)  .  In  this  latter  case,  X  is  again 
deleted,  but  we  store  x  for  future  reference. 


-  i  ■  (l  '  c,  vW ... 
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1 0 .  Boundary  points . 

The  process  just  described  in  Section  9  car.  sometimes  elimi¬ 
nate  points  of  the  boundary  of  X  0  which  lie  in  X  .  Suppose 
that  for  some  X  ,  we  find  (Xj  ,  ’’’  *xn)  >  0  for  some  i  ■ 

1,  •••  ,n  .  Then  we  can  delete  all  of  X  except  for  any  boundary 

points  of  X"  occurring  at  x-  ■  x-^  .  Any  other  boundary 

points  of  X '  which  are  in  X  are  thus  deleted. 

The  quadratic  method  of  Section  ?  deletes  any  point  x 

—  ffn 

where  fix'.  >  f  whether  x  lies  on  the  boundary  of  X  or 
not.  However,  the  Newton  method  of  Section  6  and  the  procedure 
ir.  Section  3  (which  considers  convexity)  cannot  delete  anv  bound¬ 
ary  points  of  X  '  . 

Suppose  we  applv  a  step  of  the  Newton  method  to  a  box  X 
and  obtain  a  new  box  X'  contained  in  X  .  A  simple  way  to  pro¬ 
ceed  is  to  retain  the  smallest  box  containing  both  X'  and  ail 
boundary  points  of  X  which  are  in  X  .  This  will  generally 
save  points  of  X  outside  X'  thus  reducing  the  efficiency  of 
the  procedure .  In  fact,  it  nay  be  that  the  smallest  box  containing 
the  boundary  points  of  X'p'  which  are  in  X  is  X  itself.  If 
this  were  the  case,  we  would  bypass  the  Newton  step  for  the  box  X  . 
This  approach  would  relv  upon  the  methods  of  Sections  "  and  9  to 
delete  boundary  points  of  . 

This  same  idea  can  be  used  for  the  method  of  Section  5.  If 
f  is  not  convex  in  X  ,  we  can  simply  replace  X  by  the  small¬ 
est  (perhaps  degenerate)  box,  say  T  ,  containing  the  boundary 


points  of  X  which  lie  in  X  .  In  this  case,  either  £  •  X 
or  else  X  is  a  degenerate  box  of  dimension  less  than  that  of 


X  . 

Suppose  we  are  given  a  box  X  .  For  this  approach,  if  £  ■ 
X  ,  we  do  net  applv  either  the  Newton  method  or  the  convexity 
test.  Ke  could  use  the  Newton  method  in  this  case  also  or  we 
might  bypass  its  use  whenever  X  contains  boundarv  points  of 


A  more  straightforward  procedure  is  to  simply  express  the 
boundary  of  X  as  2n  separate  (degenerate)  boxes  of  dimen¬ 
sion  n  -  1  .  The  interior  of  XlC*-’  can  then  be  treated  as  a  box 
wherein  the  global  minimum  must  be  <*  stationary  point.  However, 
the  (n-1)  -  dimensional  faces  of  X  1  have  (n-2)  -  dimensional 
boundaries  which  must,  in  turn,  be  separated  from  the  interiors, 
and  so  on.  Finallv  the  vertices  of  X'  would  have  to  be  sepa¬ 
rated.  These  vertices  alone  are  2n  points.  Even  for  moderate 
values  of  n  ,  this  separation  process  produces  too  many  (degen¬ 
erate)  boxes.  Thus  it  is  better,  in  general,  not  to  try  to 
separate  the  boundaries  from  X10"*  . 

These  two  approaches  represent  extreme  cases.  Intermediate 
methods  might  be  used  wherein  the  boundaries  of  X  in  a  given 
box  X  are  separated  off  under  special  circumstances. 

It  should,  perhaps,  be  pointed  out  that  the  Newton  method 
can  delete  boundary  points  of  X^°-  under  certain  circumstances. 


O 
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Suppose  our  algorithm  has  produced  a  degenerate  box  X  which  is 
all  or  part  of  a  face  of  X  .  In  this  degenerate  (n-1)  - 
dimensional  box,  the  Newton  method  can  delete  points  which  are 
not  in  the  (n-2)  •  dimensional  boundary  of  the  face  of  X  1  . 

Such  deleted  points  are,  of  course,  on  the  boundary  of  X  C  . 

In  some  examples,  we  shall  know  a  priori  that  the  global 
minimum  is  a  stationary  point.  In  this  case  we  are  free  tc  delete 
boundary  points  by  any  of  our  procedures. 


11.  The  list  of  boxes. 


When  we  begin  our  algorithm,  we  shall  have  a  single  box  X  Jt. 
We  appiv  the  four  procedures  described  in  Sections  5,  6,  ",  and  9 
to  this  box.  It  is  possible  that  none  of  these  procedures  can 
delete  any  of  X  '  .  If  sc,  we  divide  X 1 ^  in  half  in  a  direc¬ 

tion  of  its  maximum  width.  We  put  one  of  these  new  boxes  ir.  a 
list  L  to  be  processed  and  work  on  the  other.  These  and  subse¬ 
quent  boxes  may  also  have  to  be  subdivided,  thus  adding  tc  the 

list  L  of  boxes  ye:  to  be  processed.  If  boundary  points  are 
handled  appropriately,  the  four  procedures  described  in  Sections 
5,  b.  '  and  9  can  each  produce  more  than  one  new  sub -box;  and  all 

but  one  are  added  to  the  list.  Thus  the  number  of  boxes  ir.  the 

list  tends  tc  grow  initially. 

Eventually,  however,  the  boxes  become  snail  and  often  a  box 
is  entirely  eliminated.  Thus  the  number  of  boxes  in  the  list 
eventually  decreases  to  one,  or  just  a  few,  or  to  none  at  all  when 
c.  is  cho  on  to  be  nonzero. 
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1 2  .  Subdividing  a  box . 

In  the  initial  stages  of  our  algorithm,  we  shall  be  applying 
it  to  large  boxes.  For  example,  we  begin  by  applying  it  to  the 
entire  initial  box  Xl<^  .  Thus  it  could  be  that,  for  a  given 
box  X  ,  none  of  the  procedures  described  in  Sections  S,  6,  ", 
and  9  car.  delete  any  of  X  .  When  this  occurs,  we  wish  to  sub¬ 
divide  X  . 

We  could  subdivide  each  component  X i  of  X  into  two  parts. 

3ut  this  would  give  rise  to  2n  sub-boxes.  To  prevent  generation 

of  too  many  sub-boxes,  we  shall  divide  only  one  component  in  half. 

It  is  best  to  subdivide  the  largest  component  X^  to  prevent 

generation  of  a  long,  narrow  box. 

L  R 

Suppose  we  divide  XJ  ■  [x,  ,  x.  ]  in  half  giving  two  new 
boxes  X'  and  X’’  whose  i-th  components  are  XV  •  [xV.x.J 
and  X."  »  fx.  .x^)  ,  respecti  vely ,  where  3T  •  (xV*xjR)/ 2  . 

The  boxes  X’  and  X"  have  a  boundary  in  common  at  xi  ■  x;  . 

If  f  had  a  global  minimum  on  this  common  boundary,  we  would  sub¬ 
sequently  find  it  twice.  This  is  unlikely  to  be  the  case.  To 
avoid  having  the  same  points  in  two  boxes,  we  could  define  one  of 
them  in  terms  of  a  half  open  interval.  Thus  we  could  define 


X.  ' 
X 


[xiL,xi)  . 


It  is  simpler  to  always  use  closed  intervals.  The  extra  work 
of  keeping  track  of  whether  an  interval  contains  a  given  endpoint 
is  probably  not  worth  the  effort.  In  practice,  we  have  elected  to 
avoid  this  problem.  Thus  we  have  always  used  closed  intervals 


only.  In  general,  this  does  not  cause  the  algorithm  to  find  a 
given  global  minimum  more  than  once. 

13.  Termination . 

If  we  have  chosen  Ej  >  0  ,  we  can  continue  our  algorithm 
until  x'  '  is  entirely  eliminated.  As  pointed  out  in  Section  S, 
we  then  have  f*  bounded  to  within  a  error  t 1  .  In  this  case, 
we  do  net  obtain  a  bound  on  x*  . 

f  01 

If  Cj  •  0  ,  we  cannot  eliminate  all  of  X1  since  we  always 

retain  a  box  or  boxes  containing  the  pointful  x*  where  f'x*  « 
f*  .  As  pointed  out  above,  we  might  also  retain  a  box  or  boxes 
wherein  f  has  a  value  very  near  f*  but  no  value  equal  to  f*  . 

S.  pose  that  at  some  stage  in  our  algorithm,  the  list  L 
contains  s  boxes.  Denote  these  boxes  by  X  *  ,  •••  ,X  1  .  Let 
X,  ’  1  denote  the  interval  defining  the  i-th  component  (dimension' 
of  X<*  .  Let  w(X,  denote  the  width  of  the  interval  X. 

and  let 


w  .  [Wfx  U)l] 

i  lrisn1  J 


That  is,  w  is  the  maximum  dimension  of  X1  ■  . 

We  could  continue  processing  boxes  in  the  list  L  by  our 
algorithm  until 


3 


for  some  e,  >  0  .  This  is  provided  e,  is  chosen  large  enough 
that  the  prescribed  precision  is  attainable  using  (say)  single 
pj-^cision  arithmetic.  However,  it  is  more  convenient  computation 
ally  to  require  only  that 

(13.:)  "i  5  c2 

for  each  i  •  1,  •••  ,s  .  If  e]>  >  0  ,  we  set  c,  -  0  for  con¬ 
venience. 

Thus  whenever  a  new  box  X  ^  is  obtained  by  our  algorithm 
we  can  check  whether  (13.2)  is  satisfied.  If  sc,  we  no  longer 
apply  our  algorithm  to  X11  (except  as  discussed  below).  If 
X 1 1  contains  a  point  x*  where  f  is  a  globa*  minimum, 
then  the  location  of  x*  is  bounded.  In  fact,  if  x  is  the 

center  of  X(i)  and  (13.2)  holds  for  X(l  ,  then 

lx/11  -  x^*  5  c2/2  (:  ■  1.  •••  ,n)  • 

Let  s  denote  the  number  of  boxes  remaining  and  denote 
the  boxes  by  X^11  (i  -  1,  s) .  As  a  final  step,  we  want 

to  assure  that  f *  is  bounded  sufficiently  sharply,  *e  do  this 
as  follows. 


o 
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For  *3Ch  i  ■  1,  *'•  ,s  we  evaluate  f(X  *  ’)  ;  that  is  we 
evaluate  f  with  interval  arguments  X^i'*  (j«l,  •'*  ,n)  . 

The  result,  say  (F , L , F . R]  contains  the  range  of  f  for  all 
x  <  X^  ,  but  will  not  be  sharp,  in  general.  However,  if  t, 
was  chosen  tc  be  snail,  the  interval  result  should  be  ’'close  to 
sharp"  since  c,  is  an  upper  bound  on  the  largest  dimension  of 
any  box  X(l)  ;  and  the  smaller  X  1  is,  the  sharper 
[F.L,F.R]  is  .  (See  [9].'  Therefore,  it  is  generally  not 
necessarv  to  use  special  procedures  tc  sharpen  the  computed  inter¬ 
val  . 

Since  contains  the  range  of  f(x)  for  all 

x<r  ,  we  have 

Fj1  5  f(x)  s  FiR 

for  anv  x  <  X* 1  .  Denote 

,  ruin  r  L 

-  •  15155  Fi 

Then 

f  5  f(x) 

for  any  x  in  any  of  the  boxes  X '  *  (i  ■  1 ,  ’  *  s )  •  iheretcre, 

since  any  global  minimum  must  occur  at  a  point  x*  lying  in  one 
of  the  boxes  ,  we  have  f  5  f*  .  But  also  f*  £  I  (as 

discussed  above)  and  hence 

)  f  s  f*  ;  f  . 
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We  thus  have  bounds  on  f*  .  However,  they  may  not  be  sharp 
enough  since  our  specified  requirement  is  to  bound  £*  to  within,  sa 


e 0  ;  and  it  may  be  that  f 


f  >  e, 


If  this  is  the  case,  we 


shall  improve  our  bounds.  To  do  this,  we  find  a  value  i  for 
which  f  »  f.“  .  v.'e  applv  our  main  algorithm  to  X  •  .  This 

will  increase  £  ^  ,  in  general.  It  might  also  decrease  T  . 

For  exact  interval  arithmetic,  this  must  decrease  f  •  f,“  since 
X  -  will  be  reduced  in  sice  '’ever,  if  it  is  merely  subdivided' . 
Repeating  this  step  for  each  i  such  that  f  »  £.~  ,  we  must 
decrease  f  •  f  (at  least,  if  exact  interval  arithmetic  is  used 
and  hence  eventual lv  have 


f  fi  e. 


so  that  f*  is  bounded  to  sufficient  accuracv  since  (13.3'!  holds. 

Because  of  rounding  errors,  we  cannot  reduce  T  -  f 
arbitrarily,  in  practice.  Hence  we  assume  is  chosen 
commensurate  with  achievable  accuracy  using  (say)  single 
precision  arithmetic. 

He  also  require  that 

»i*  •  ^  1  ‘3 

‘for  each  box  X;i>  (i  •  l,***,  s).  For  convenience,  we  can 
choose  t .  ■  e0  to  reduce  the  number  of  quantities  to  be  specified. 
Note  that  F^  <_  7  since  otherwise  f(x)  >  7  for  all  x  i  X1  i 
in  which  case  Xu  can  be  deleted.  Hence 

FiR  I  *  ♦  c3  <  f  *  cQ  *  c3  <  f*  ♦  c0  *  c3 
for  every  i  •  l,***,  s.  That  is,  every  remaining  box  contains 
a  point  x  at  which  f(x)  differs  from  f*  by  no  mere  than 


1 4  .  The  steps  of  the  algorithm . 

He  now  describe  the  steps  involved  in  our  algorithm.  Initially 
the  list  L  of  boxes  to  be  processed  consists  of  a  single  box 
X1'1'1  .  In  general,  divide  the  list  l  into  the  list  L ^  of 
intervals  X1,  satisfying  the  condition  s  c,  (see  (13. 1; 
and  a  list  L,  which  do  net  satisfy  this  condition. 

(01 

He  assure  we  have  evaluated  f  3t  the  center  of  X ‘  J  and 
thus  obtained  ar.  initial  value  for  f  .  The  subsequent  steps  are 
to  be  done  in  the  following  order  except  as  indicated  by  branching 

(1)  Of  the  boxes  in  L,  ,  choose  one  which  has  been  in  L, 
longest.  Call  it  X  .  If  L,  is  empty,  go  to  step  (11)  if 

e,  >  0  .  If  L,  is  empty  and  e .  •  0  choose  a  box  which  has  bee 
in  L,  longest  and  go  to  step  12.  If  both  Li  and  L,  are 
empty,  print  the  bounds  f  -  e.  and  f  on  f*  and  stop. 

(2)  Check  for  monoton i c l ty .  Evaluate  g(X)  as  described  in 

Section  9.  For  i»l,***,n,if  g 1 ( X )  >  0  (  <  0)  and  the 

boundary  of  X  at  )  does  not  contain  a  boundary 

point  of  X!C  ,  delete  X  and  go  to  step  (1  .  Otherwise,  if 

gx(X)  >0  (50)  ,  replace  Xt  -  [x^.x^J  by  [x^.x,1] 

R  R 

((Xj  .x/j)  .  Rename  the  result  X  again. 

(3)  Test  for  non  -  convex! ty  as  in  Section  5.  Let  X’  denote  the 

smallest  box  in  X  containing  all  the  boundary  points  of  X 
which  lie  in  X  .  If  X’  •  X  go  to  step  4.  Otherwise,  evaluate 
h^fX,,  **„)  f°r  *  "  *’*  •  If  the  resulting  interval 

is  strictly  negative  for  any  value  of  i  ,  replace  X  by  X’ 
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If  X’  is  empty,  go  to  step  1  .  If  X’  is  not  empty,  put  it 
in  the  list  L  and  go  to  step  1. 

(4)  Begin  use  of  the  quadratic  method  of  Section  For  those 
values  of  i  ■  1,  •••  ,n  for  which  the  left  endpoint  of  H^(X) 
is  non-negative,  solve  the  quadratic  for  the  interval  Y^  to 
replace  X^  .  Rename  the  result  Xj  . 

(5  Begin  use  of  the  Newton  method.  For  those  values  of  i  • 

1,  •••  ,n  for  which  0  i  [BJ  X^]^  solve  for  the  new  interval  to 
replace  X^  .  Rename  the  result  Xj  .  For  a  given  value  of  i  , 

omit  this  step  if  a  reduction  of  X.  will  delete  boundarv  points 
cf  X'  from  the  box  X  . 

(6)  Complete  the  quadratic  method.  For  these  values  of  i  not 
used  in  step  4,  solve  the  quadratic  for  Y,  .  If  Yx  is  a  single 
interval,  replace  X,  bv  Yj  ,  renaming  it  X.  .  Otherwise  save 
Yi  for  use  in  step  8. 

(*  Complete  the  Newton  method.  For  those  values  of  i  not  used 
in  step  3,  solve  for  the  new  set  (say)  YV  .  For  a  given  value 

of  i  ,  omit  this  step  if  a  reduction  of  Xj  will  delete  bound¬ 

ary  points  of  X'u  from  the  box  X  .  If  YV  is  a  single  inter¬ 
val,  replace  by  YV  ,  renaming  it  Xi  . 

(8)  Combine  the  results  from  the  quadratic  and  Newton  methods 
for  those  components  Xi  for  which  both  methods  divided  X^  into 
two  sub- intervals .  That  is,  find  the  intersection  Yj"  of  Y. 
from  step  6  and  Y^ ’  from  step  ?.  If  Yi"  is  composed  of  three 
intervals,  replace  it  by  either  Y ^  or  Y^’  ,  whichever  has  the 


3" 


smallest  intersection  with  Xj  .  Of  all  the  Y^"  ,  save  the  one 
(or  two  or  three'  which  deletes  the  largest  subinterval  of  X,  . 
That  is,  save  that  Y  "  whose  complement  in  Xj  is  largest. 

Let  i  be  its  index.  For  all  relevant  values  of  i  *  i  ,  re¬ 
place  Y/'  by  X (  ,  that  is,  ignore  the  fact  that  part  of  X; 
could  be  deleted. 

(9)  If  Y  ”  exists.  That  is,  if  at  least  one  interval  was 

divided  into  two  sub- interval s ,  say  Y.  and  Y.  ,  subdivide 

»  J 

the  box  X  into  two  sub-boxes.  These  sub-boxes  will  have  the 
same  components  Xi  as  X  except  one  will  have  j-th  component 
Y ,  *  and  the  other  will  have  j-th  component  Y . 1  *  .  If  no  such 

Y  "  exists,  we  nav  wish  to  subdivide  the  current  box.  * 

denote  the  box  chosen  in  step  1  and  let  X”  denote  the  current 

box  resulting  from  applying  steps  2  through  8  to  X  .  * f  the 

improvement  of  X"  over  X  is  so  small  that  (sav) 

v(X")  >  . “5w(X) , 

then  divide  X”  in  half  in  its  greatest  dimension. 

,1C  Evaluate  f  at  the  center  of  the  box  or  boxes  resulting 
after  step  9.  Update  T  as  described  in  Section  4.  Put  the 

box ( es)  in  the  list  L  and  go  to  step  1. 

(11)  Evaluate  f(X'^)  for  each  remaining  box  X  1  in  L. 

‘‘Denote  the  result  by  [F^.F^].  If  F^  -  F^  >  cQ  for  anv 

value  of  i,  use  Xu  for  X  and  go  to  step  4.  If  Fi  •  Fj  i 

for  all  i  ■  s,  f ind 

f  .  miTi  f  1 
1  1liLs  1 

Then  print  the  bounds  f  and  T  on  f*  and  stop. 


For  c,  >  0,  these  steps  bound  i*  to  within  an  error  e 


•  0.  thev  found  f*  to  within  eA;  tnev  bound  x*  to  within 
1  ’  o 

and  thev  assure  that  for  any  point  x  in  any  final  box,  f(x' 


exceeds  f*  bv  no  more  than  •: 


es  branch  to  step  i 


we  seme 


we  could  go  to  step  2,  but  it  is  unlikely  that  either  step 
step  3  will  be  helpful.  This  is  because  we  expect  each  box 


re’-ainmg  at  this  stage  to  contain  a  minimum. 


15.  A  numerical  exa 


he  now  illustrate  the  steps  of  our  algorithm,  he  shall  cor. 
sider  the  so-called  three  hump  camel  function. 


has  components 


The  interval  Jacobian  J(X)  (see  Section  3)  has  elements 


As  described  in  [4],  a  better  formulation  for  J(X)  coulc 
be  derived  which  would  give  smaller  intervals,  in  general.  How¬ 
ever,  we  shall  use  the  simpler  form  given  here. 

Suppose  that  the  box  we  choose  in  step  1  has  first  component 
v  -  11  l.ll.  In  step  2  we  find  that,  whatever  X,  is, 


Since  this  is  strictly  negative,  we 


Hence  i 


for  anv  value  of  x 


minimum  in  the  interval  X 


does  not  contain  3  boundary  point  0 


the  box  chosen  in  step  1  has  compcner. t s 


ve ,  we  cannot  sav 


convex  in 


we  evaluate 


is  non-positive  for  all  x*X  and  hence  f 
Thus  we  can  replace  X  by  the 


We  see  that 

is  smallest  in  X  for  x2  -  1 
degenerate  box  X'  with  components  X 
If  the  box  X  had  had  components 
12.31  .  we  would  have  obtained 


[I] 


In  this  case  g,(x)  is  strictly  positive  and  we  can  eliminate  all 

to 

of  X  unless  the  boundarv  of  X  at  x,  ■  2  contains  a  boundarv 
point  of  .  "Suppose  v  •  [0,1]  and  X-/0"1  ■  [1,3]  . 

Then  X  ^  has  boundary  points  at  x.  *  2  for  x1  •  0  and  1  . 
We  could  thus  delete  all  of  X  except  the  points  (0,2)  and 


(1,2)  .  This  is  simple  to  do  in  this  two-dimensional  problem. 

In  higher  dimensions,  it  might  be  simpler  to  retain  the  entire 
boundary  at  x,  •  2  . 

Now  suppose  that  X  is  given  bv  X^  *  X,  ■  [0,1]  .  Then 
g,(X)  ■  [-5.2,5]  and  g,(X)  -  [-1,2]  so  that  we  dc  not  have 
monotonic: tv .  Therefore,  we  dc  step  4  which  involves  the  quadra 
tic  method  of  Section 

For  this  box,  we  obtain 


Hji (X)  •  J1:(X)  -  [-8.6,9]  .  H;1(X)  •  2J;1(X)  -  -2  ,  H;;(X)  »  J;;(X^  -  2 


The  center  of  the  box  is  at  x  *  ( . 5 ,  .5  .  We  wish  to  evaluate 
f (x)  and  g(x)  .  W'e  cannot  obtain  f(x)  exactly  using  finite 
precision  decimal  arithmetic.  Let  us  use  five  significant  decimal 
digits  and  evaluate  f(x)  using  interval  arithmetic  to  bound 
rounding  errors.  Thus  we  replace  the  coefficient  1/6  by 
[.16666,  .1666"]  and  obtain 


f(x)  -  [.4369",  .43699]  . 

We  also  obtain 


g(x) 


[1.0062,  1.0063]  ' 
.5 


Suppose  we  have  previously  obtained  f 


and  that  we 


choose  *  0  .  To  do  step  4,  we  wish  to  solve  (7.2)  for  points 

v  <  X  where  we  know  that  f(y)  >  f  does  not  hold  and  hence 
f(y)  5  7  night  hold.  If  we  first  tried  to  solve  for  Yj  ,  we 
would  rewrite  (7.2)  in  the  form  (7.4).  However,  the  left  endpoint 
of  is  less  than  :ero.  Hence,  solving  C-4)  would  give 

rise  to  two  semi - inf inite  intervals  Therefor*,  w*  defer  this 
operation  until  step  6  and  first  solve  for  Yj  which  will  be  a 
single  interval  since  the  "left  endpoint"  of  H^C*)  positi\e 

We  solve  for  Y,  using  C.5).  As  we  have  not  yet  solved  for 
Yj  ,  we  use  X1  in  its  place.  Substituting  into  (7.5),  we 
obtain 

(-  1.2412  ,  1.9652)  *  ( 0 , 1  ]  >' ,  *  V;  i  0. 

From  equation  (7.8),  the  solution  set  is 
2  2  ■  ( -1  -  "21 2  ,  1-1141)  . 

Hence 

2,  •  x,  ♦  Z,  •  (-1.2212,  1.6141) 

and  ^2  "  “2  '  *2  "  *2  ’ 

Thus  we  have  not  deleted  any  of  X,  • 

Next  we  do  step  5  which  app lies  that  part  of  the  Newton 

‘method  which  generates  a  single  interval.  The  interval  Jacobian 

f  ( * 8 . 6 ,  9)  -1*1 


The  center  of  this  interval  matrix  is 


;  .  P  ^ 

c  l-i  :J 


whose  inverse  is  (approximately) 


f  -3.3333  -1.666'  1 


L  -1.666" 


.  oocoo  j 


Fcr  simplicity  of  exposition,  we  shall  compute  EJ(X) 
explicitly,  he  obtain  from  (6.2)  (with  ;  replaced  by  X  )  , 

f[-4. 18",  -4.18-211  r  [-28.334,  30.334]  -.0001  1 

(13.4)  1  ♦  Cy-x)  ■  0  . 

[1-1.844,-1.8436]  J  J-14. 668,  14.663]  [1,  1.0001]  j 

he  try  to  improve  X.  first  rather  than  X,  because  [BJ  X  1^, 

contains  :ero  while  [BJ(X)],,  does  not.  Thus  the  first  equation 

*  * 

of  (13.4)  gives  rise  to  two  new  intervals  while  the  second  equa¬ 
tion  does  not. 

The  second  equation  is 

[-1.844,  -1.S436]  *  [-14.668,  14.668]  (X1-x1)  *  [1 ,  1.0001]  (y2*x:)  ■  0 

-.where  we  have  replaced  y^  by  Xj  .  Solving  for  .  we  obtain 


the  interval 


[-4.9904,8.6"8]  . 


O 


The  intersection  Y,  »  2,  h  X,  equals  X?  so  again  no  improvement 
has  been  made . 

Step  o  prescribes  that  we  use  the  quadratic  method  to  try  to 
improve  those  components  of  X  not  solved  for  in  step  4.  Ke 
solve  (".S)  for  points  where  f(y)  5  i  .  We  would  use  Y, 

in  place  of  X,  ;  but  they  are  equal.  Substituting  into  (”.5), 
we  obtain 

[.0869",  . 8  3c  99 ]  ♦  (.  506:.  1.5063]y\  *  [  -  4 . 5 , 4 . 5  ]  y  x  2  <  0 

Using  equation  (".11),  we  find  that  this  quadratic  has  the 
solution  set 

2 ,  •  ['*,  -  .10093]  c  [ . 4  2  S  5  5  .  ®  ]  . 

Thus 

:1  "  "I  *  X1  '  l*“*  -3990"]  J  [.92555,  «] 

and 

Y1  "  :1  "  X1  "  -3990"]  -  (.92555,  1] 

note  we  have  eliminated  a  subinterval  of  length  .52o4S  from  Xj  . 

In  step  ",  we  use  the  Newton  method  to  try  to  improve  X^  . 

We  solve  the  first  equation  of  (15.4), 

[-4. IS"",  -4.18”2]  *  [-  28.334  ,  30.334  ]  (Vj-Xj)  -  .0001  (Xj-x,)  ■  0, 

where  we  have  now  replaced  Y,  by  X,  •  Solving  for  Vj  ,  we 
obtain  the  two  semi  -  inf inite  intervals 

:1(3)  ■  [  — ,  .  35223],  21(4)  -  [.63803.-]  • 

Their  intersections  with  X1  are  (say)  Y^31  and  Y^4^  where 

Y^3^  -  [0, .35223]  ,  Yj(4)  • 


[.63803,1]  . 
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Ke  wish  to  combine  the  results  obtained  using  the  quadratic 
method  and  the  Newton  method.  Thus  we  retain  the  intersection  of 

Yj(1)  J  Y^21  and  Y^3^  U  ■ 4 1  which  is 

[0,  . 35223 )  j  (.  92555  ,  1  ]  . 

Ke  have  deleted  a  substantial  portion  of  the  original  box  X  . 

The  remaining  points  compose  the  two  boxes 


"(0,. 

35223)' 

'(.92555  ,  1  y 

-  10. 

1)  - 

and 

[0,1] 

In  subsequent  steps,  our  method  would  be  applied  separately  to 
each  of  these  boxes  . 

1 b .  Computational  results 

We  now  describe  some  cor.putat  ior.al  results  obtained  using  the 
algorithm  described  above.  The  computations  were  done  or.  the 
Amdahl  4“0V/6-II  computer.  In  each  case,  we  assumed  it  was  known 
that  the  global  minimum  occurred  in  the  interior  of  the  initial 
box.  This  speeds  up  the  algorithm  since  boundary  points  need  not 
get  special  treatment. 

This  is  consistent  with  the  fact  that  we  are  really  considering 
the  unconstrained  case.  Ke  intend  to  treat  the  constrained  case 
in  a  later  paper. 

V«  give  results  for  only  one  example  which  typifies  the  problems  in  two  and 
three  dimensions  which  we  have  used.  The  example  Is  the  three  hump  caael  function 
given  by  (.15.1). 

This  function  has  its  global  slalma  at  the  origin.  It  has  two  local  minima 
at  approximately  1.75,  ♦  0.8^  and  two  saddle  points  at  approximately  1.07, 

±  0.533).  Out  initial  box  was  deflnad  by  •  Q-2,  Zj  which  contains  all 


these  points.  We  chose  f ^  •  0  and  t ^  •  f ^  •  10“*. 


- 


We  find  that  after  eight  steps  of  our.  algorithm,  we  have  six 
sub-boxes  in  our  list.  In  the  next  step  a  sub-box  is  entirely 
eliminated  and  after  the  fifteenth  step,  only  one  sub-box  remains 
but  its  width  exceeds  C2»  After  an  additional  step,  we  obtain 
the  final  box 

’[-5.78  x  10'7,  7.11  x  10'V 

X  *  [-2.91  x  10*7,  3.56  x  10'6] 

L  J 

Here  and  in  the  following,  we  record  results  to  only  three 

decimals.  This  box  satisfies  the  error  criterion  requiring 

its  width  to  be  less  than  c?.  Evaluating  f  at  the  center  of 

this  box,  we  obtain  f  ■  1.12  x  10 

As  prescribed  in  step  ll,  we  evaluate  f(X)  and  obtain 

[-1.24  x  10*11,  7.57  x  1 0 ’ 1 °] .  Thus  f  --1.24  x  10*11  and 

f*  (  [f  ,  f]  -  [-1-24  x  10*U,  1.12  x  10*10]. 

Since  f  -  f  <  cQ,  we  have  f*  bounded  to  the  prescribed 
tolerance.  If  we  approximate  f*  by 

(f  ♦  f ) / 2  ■  4.98  x  10"11, 

then  we  know  that  the  error  is  at  most  6.22  xlO  11  in  magnitude. 

If  we  approximate  x*  by  the  center  (3.27  x  10  1.63  x  10 

then  we  know  that  the  error  in  Xj*  is  less  than  3.85  x  10  6  and 
the  error  in  X2*  is  less  than  1.93  x  10 

We  have  obtained  x*  to  far  more  accuracy  than  required 
because  of  the  rapid  rate  of  convergence  of  the  interval  Newton 
method  used.  The  bound  on  f*  is  much  better  than  required 
simply  because  a  given  error  bound  on  x*  automatically  yields 
a  much  better  bound  on  f*  . 


We  also  used  this  example  with  an  initial  box  of  width  2  x  IQ6.  This  case 
required  46  steps  to  run  to  completion.  This  illustrates  that  If  we  use  a  very 
large  box  to  assure  containment  of  x* ,  the  computing  time  need  not  increase 
drastically. 

1  “  .  Corel  us  icr.  ■ 

he  have  presented  an  algorithm  for  solving  the  unconstrained 
minimisation  probier  assuring  we  have  an  initial  box  which  is 
known  tc  contain  the  minimum. 

It  would  certainly  be  possible  to  construct  a  highly  oscillato 
function  for  which  our  algorithm  would  be  prohibitively  slow. 
However,  it  has  converged  adequately  rapidly  for  the  test  problems 
or.  which  we  have  tried  it.  (See  Section  16.) 

ke  have  assumed  f(.x)  t  C “ .  The  global  nir.imicat  ion  problem 
car.  also  be  solved  for  f:..\  (  C1 .  In  this  case,  the  Newton 
method  cannot  be  used.  The  quadratic  method  car.  be  replaced  by  a 
corresponding  linear  method  in  which  we  find  points  y  at  which 
f (y)  <_  7.  This  is  done  by  noting  that  if  x  t  X  and  y  <  X,  then 

f (y)  -  f(x)  ♦  (y  -  x)rg(0 

for  some  £  <  X.  Thus  we  can  solve  for  the  approximate  point?  v 
from 

f(x)  *  (y  -  x)T  g(X)  <  7 

For  the  problems  of  low  dimension  on  which  we  have  used  this  metho 
it  was  less  efficient  then  the  quadratic  method  described  in 
Section  7.  Ke  do  not  know  the  relative  efficiencies  for  large  n. 


It  is  possible  to  solve  the  global  optimitation  case 
when  f'x)  is  only  continuous  but  not  differentiable.  However, 
our  algorithm  is  very  slow.  It  entails  a  different  approach 
that  we  hope  tc  describe  in  another  paper. 

The  nonlinear  constrained  optimization  problem  car.  also 
be  solved  by  interval  methods.  An  extension  of  our  algorithm 
is  required.  Our  experience  in  this  case  is  for  hand  caiculat 
only.  A  difficulty  exists  (currently)  when  it  is  difficult  to 
find  a  point  in  the  neighborhood  of  x*  which  is  without  quest: 
feasible . 

One  of  the  virtues  of  interval  arithmetic  is  that  it  is 
usually  possible  to  formulate  an  iterative  algorithm  in  such  a 
way  that  it  stops  automatically  when  the  best  possible  result 
has  been  obtained  for  the  finite  precision  arithmetic  used. 

Ke  plan  to  do  this  for  our  algorithm  and  thus  preclude  the 
need  for  specifying  t ^ ,  c ^ ,  c,,  and  e.. 


o 
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